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Final  Report 

Application  of  the  Oriented-Eddy  Collision 
Model  to  Complex  Turbulent  Flows 
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Contract  Number:  N0001 4-08-1  -0275 
Dates:  1/1/08-12/31/10 
Program  Officer:  Ronald  Joslin 

Graduate  Students:  Michael  B  Martell  Jr 

Summary 

The  Oriented-Eddy  Collision  (OEC)  turbulence  model  hypothesizes  that 
turbulent  flow  can  be  modeled  as  a  collection  of  interacting  fluid  particles  (or  eddies) 
which  have  fluctuations  and  an  inherent  orientation.  The  model  has  been  formulated  in 
the  form  of  a  set  of  partial  differential  equations.  Underlying  this  approach  is  a 
unique  PDF  collision  model  that  includes  the  effect  of  orientation  information  along 
with  the  usual  position  and  velocity  information  in  the  formulation  of  the 
probability  density  function.  This  adds  important  physics  to  the  model  and 
differentiates  it  from  most  other  PDF  models  and  Reynolds-Averaged  Navier- 
Stokes  models. 

The  Oriented  Eddy  Collision  model  exactly  captures  rapid  distortion,  which 
is  a  major  shortcoming  of  most  prior  Reynolds  stress  transport  models.  The  ability 
to  predict  highly  non-equilibrium  flow  situations  well  is  a  major  feature  of  the  model. 
The  model  automatically  (via  its  construction)  satisfies  realizability  and  other  known 
mathematical  constraints.  It  is  readily  extensible  to  complex  geometries  and  additional 
physical  affects  (such  as  compressibility,  particles,  etc).  The  model  has  been 
implemented  in  the  open  source  CFD  framework  OpenFOAM  for  rapid  dissemination. 
It  has  been  tested  extensively  on  basic  turbulent  flow  benchmarks  and,  more  recently, 
on  flows  with  solid  boundaries. 

Background 

The  traditional  approach  to  modeling  turbulence  (or  other  types  of  non- 
Newtonian  fluids)  is  to  hypothesize  equations  for  the  unknown  stress  tensor,  lin 
turbulence  this  is  the  Reynolds  stress  tensor.  Due  to  the  fact  that  the  eddies  which 
make  up  the  flow  are  roughly  the  same  size  as  the  gradients  in  the  mean  flow  these 
eddies  respond  on  similar  timescales  as  the  mean  flow.  This  means  that  algebraic 
models  are  rarely  predictive,  and  time-dependent  evolution  equations  for  the  stress 
tensor  must  be  hypothesized.  In  turbulence,  these  evolution  equations  are  the  exact 
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but  unclosed  Reynolds  stress  transport  (RST)  equations.  Simpler  turbulence  models, 
such  as  the  k-e  model  or  algebraic  Reynolds  stress  models  are  simplifications  of  the 
RST  equations. 

There  is  a  strong  analogy  between  turbulent  fluid  flow  and  non-Newtonian  or 
granular  flows.  Very  similar  to  turbulent  flows,  transport  equations  are  very  often 
developed  for  non-Newtonian  stress  tensors  (the  Oldroyd-B  model  and  FENE-P 
models  (Herrechen,  et  al.  1997)  are  examples).  In  fact,  we  note  that  many  important 
turbulence  modeling  concepts  (realizability,  material  frame  indifference,  tensor 
consistency)  actually  find  their  origins  in  the  non-Newtonian  literature  at  this  transport 
equation  level.  This  work  is  predicated  on  treating  turbulence  modeling  in  a  fashion 
that  is  similar  to  non-Newtonian  fluid  modeling. 

It  has  long  been  recognized  in  the  non-Newtonian  fluid  community  that  transport 
equation  models  have  serious  limitations.  An  alternative  approach  is  to  model  the  fluid 
at  the  particle  collision  level  rather  than  using  a  transport  equation  for  the  stress.  This 
approach  is  more  versatile,  and  in  many  ways,  more  fundamental.  For  example, 
modeling  a  gas  as  particles  with  binary  elastic  hard  sphere  collisions  gives  the  Navier- 
Stokes  equations  and  the  perfect  gas  law  when  the  density  is  high,  but  also  the  correct 
gas  behavior  even  when  the  density  is  low  (when  Navier-Stokes  is  not  valid).  In  this 
work,  we  investigated  the  possibility  of  modeling  turbulence  as  a  collection  of 
interacting  oriented  particles  (which  will  turn  out  to  be  disks  or  rods). 

Once  a  certain  collision  behavior  has  been  hypothesized  there  are  three  very 
different  ways  to  solve  the  particle  system  numerically  and  obtain  a  prediction  of  the 
fluid  behavior.  The  most  straightforward  technique  is  the  ‘molecular  dynamics’ 
approach  where  one  numerically  tracks  all  the  particles  in  the  domain,  and  performs 
collisions  when  they  occur.  This  approach  has  a  computational  cost  equivalent  to  large 
eddy  simulation  (LES)  and  is  not  considered  further.  The  other  two  approaches  note 
that  one  does  not  really  care  what  happens  to  individual  particles  but  only  what 
happens  to  particles  on  average.  The  quantity  of  interest  then  becomes  the  probability 
density  function  that  describes  the  probability  that  a  particle  (at  a  certain  place  and 
time)  has  a  certain  velocity.  The  evolution  of  the  probability  distribution  function,  f, 
obeys  the  exact  equation 

df  df  df  dn  df  df 

—+v,—+a.—+ — L— =  —  m 

dt  dx.  dvt  dt  dnt  dt  collisions 

where  a,  is  the  acceleration  due  to  external  forces  (like  gravity),  n,  is  the  particle 

orientation,  and  the  right-hand  side  describes  the  average  affect  of  the  collisions  on  the 
PDF.  It  is  this  average  collision  behavior  that  we  now  wish  the  model  to  predict.  This 
collision  term  is  also  what  gives  the  ensuing  model  its  name.  Our  collision  models 
assume  the  collision  term  has  a  Fokker-Planck  form  (see  Equations  (2)  through  (4)). 

There  are  three  different  ways  to  solve  this  PDF  equation.  Using  the 
equivalence  between  the  Fokker-Planck  equation  and  the  Langevin  equation 
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(Brownian  motion),  it  is  possible  to  construct  a  Lagrangian  particle  method.  This  is 
essentially  the  approach  extensively  researched  by  Pope  (1994,  2000)  and  coworkers 
(with  the  major  difference  from  this  work  being  that  that  Pope  and  others  typically  do 
not  use  oriented  eddies,  just  colliding  spheres).  In  this  particle  approach,  the 
Lagrangian  particles  move  like  Brownian  dust  particles.  They  move  with  the  mean  flow 
and  are  randomly  perturbed  using  a  prescription  given  by  the  model.  In  this  way  each 
particle  is  independent  from  all  the  others,  and  simply  interacts  with  the  average  of  all 
the  other  particles  (/.e.  the  mean  flow,  and  average  turbulence  statistics).  This  Monte- 
Carlo  approach  is  less  expensive  than  tracking  and  implementing  individual  collisions 
(‘molecular  dynamics’  approach)  but  is  still  expensive  because  a  large  statistical 
sample  of  particles  that  is  required.  These  methods  significantly  over-resolve  the 
shape  compared  to  what  is  necessary  to  model  the  turbulence. 

Using  a  mesh  based  method  (rather  than  Monte-Carlo  sampling)  is  also  possible.  A 
mesh  based  method  can  use  a  very  coarse  mesh,  thereby  lowering  the  costs  involved. 
A  very  coarse  mesh  in  velocity  space  is  an  idea  borrowed  from  Lattice-Boltzmann 
numerical  methods  for  solving  the  Navier-Stokes  equations.  These  methods  solve  a 
PDF  equation  with  a  very  simple  collision  term  that  is  intended  to  give  Navier-Stokes 
(Newtonian)  fluid  behavior.  The  difference  in  this  work  is  that  we  solve  a  PDF 
equation  with  a  much  more  complex  collision  term  (Fokker-Plank),  which  results  in 
RANS  behavior  for  the  fluid  (rather 
than  Newtonian).  The  coarse 
mesh  is  acceptable  in  both  cases 
because  the  interest  is  not  in  the 
PDF  itself  but  in  its  lowest  order 
moments  -  the  mean  flow  and  the 
stresses.  These  low-order 
moments  can  be  reasonably 
extracted  from  a  very  coarse 
approximation  of  the  PDF. 

Note  that  the  Langevin  approach  is 
equivalent  to  approximating  the 
PDF  with  a  random  sample,  and  a 
large  sample  is  needed  even  to 
approximate  the  low  order 
moments  reasonably  well.  The  Langevin  approach  is  slower  because  it  provides  more 
information  (about  the  higher  order  moments).  Unfortunately,  one  has  little  interest,  in 
engineering  turbulence  models,  in  the  extra  information  the  Langevin  solution  method 
provides.  These  first  two  approaches  to  solving  the  collision  model  (and  the  brute 
force  LES-type  approach  of  tracking  actual  eddies),  is  shown  in  Figure  1. 

While  the  coarse  mesh  approach  is  inspired  by  the  success  of  lattice-Boltzmann 
numerical  methods,  the  approach  must  be  numerically  different.  This  is  because  the 
PDF  governing  molecular  interactions  (Lattice-Boltzmann)  has  a  variance  (width)  that 
is  much  larger  than  the  mean  and  which  is  essentially  constant  (related  to  the  speed  of 
sound).  In  contrast,  the  PDF  for  turbulence  has  a  variance  which  is  much  smaller  than 


Collision  Models 


PDF  Methods 


Particle 

Tracking 

(‘molecular 

dynamics') 


Coarse  Discretization  Langevin  Equation 
(“lattice  methods')  (Lagrangian  particle  methods) 


Figure  1:  Taxonomy  of  classic  collision  model  approaches. 
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the  mean  (turbulence  intensities  are  measured  in  percent),  and  which  can  vary 
significantly  (in  time  or  space).  This  is  illustrated  in  Figure  2: 


Figure  2:  Left  -  a  typical  PDF  for  molecules.  Right  -  A  typical  PDF  for  turbulence. 

To  capture  the  turbulence  PDF  with  only  three  points  it  is  necessary  to  have  a  moving 
adaptive  mesh  in  velocity  space.  In  order  to  avoid  losses  due  to  interpolating  one 
mesh  to  another  as  the  mesh  moves,  we  implemented  a  fully  conservative  scheme  in 
which  the  mesh  moves  continuously  in  time  (during  the  time  step).  This  uses 
technology  previously  developed  by  Perot  &  Nallapati  (2003)  for  moving  meshes  in 
physical  space.  In  actual  practice  the  PDF  is  three-dimensional.  An  isosurface  of  an 
actual  PDF  (the  50%  value)  is  shown  below.  This  PDF  is  modeling  the  behavior  of  the 
Le  Penven  et  al  (1985)  return-to-isotropy  Case  III  >  0  experiment.  Note  the  fairly  large 
changes  in  the  shape  and  size  of  the  distribution  even  for  this  simple  experiment.  It 
can  also  be  seen  in  this  figure  that  a  spherical  PDF  corresponds  to  isotropic 
turbulence. 


Figure  3:  Evolution  of  the  50%  isosurface  of  the  PDF  for  the  return-to-isotropy  experiment  of 
Le  Penven,  et  al. 


The  coarse  mesh  approach  was  used  in  our  initial  ONR  work,  but  the  current 
model  actually  uses  a  third  PDF  equation  solution  approach  that  ends  up  being  far 
more  familiar  to  the  CFD  users.  In  the  current  project  we  took  moments  of  the  PDF 
collision  equation  over  velocity  space  in  order  to  construct  a  set  of  partial  differential 
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equations  for  the  behavior  of  each  type  of  oriented  disk.  No  velocity  space  meshing  (or 
Monte-Carlo  sampling)  is  now  necessary.  This  third  approach  can  be  directly  included 
into  PDE  code  frameworks  (such  as  OpenFoam).  It  has  the  disadvantage  that  the 
third-order  correlations  must  now  be  modeled  (in  the  pure  PDF  solutions  they  can  be 
deduced  from  the  PDF).  This  means  the  turbulent  diffusion  processes  must  now  be 
modeled. 

Theoretical  background  for  the  un-orientated  eddy  collision  (EC)  model 

Lundgren  (1967)  first  derived  the  exact  expression  for  the  collision  term  in  the 
PDF  evolution  equation  for  turbulence.  As  might  be  expected,  this  collision  term 
cannot  be  expressed  solely  in  terms  of  the  PDF,  and  solution  of  the  PDF  evolution 
equation  therefore  requires  a  model  for  the  collision  term.  Original  development  of  the 
OEC  model  focused  on  generalizations  of  the  Fokker-Plank  collision  model  (that  was 
derived  to  describe  Brownian  motion).  In  its  simplest  form  this  collision  model  has  the 
form, 


df 


dt 


collision 


(2) 


where  u^^vjdv  is  the  mean  velocity  and  a  and  b  are  model  constants.  For 

turbulence  this  needs  to  be  generalized.  Pope  and  coworkers  (Pope  2000,  Reynolds 
1995,  Van  Slooten  1996)  use  the  form, 


£ 

dt 


collision 


dV; 


\Gv'  f]+b^+vil 

L  ,y  7  J  dv,2  dx2 


(3) 


where  vV  =  v,  -  w,  is  the  fluctuating  velocity  and  the  first  term  (the  drift  term)  now  has  a 

matrix  model  parameter  Gy ,  and  a  viscous  term  has  been  added  for  near  wall  (low  Re 

number)  calculations.  The  conversion  of  these  Fokker-Planck  models  to  a  Langevin 
equation  for  numerical  solution  dictates  that  the  diffusion  term  (with  b)  be  isotropic  and 
not  have  a  tensor  coefficient. 


Original  development  of  the  OEC  model  analyzed  the  following  even  more 
generalized  Fokker  Plank  model. 
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The  last  term  on  the  right  hand  side  accounts  (exactly)  for  the  coarse  mesh  motion  in 
velocity  space  (to  track  the  PDF).  The  first  three  terms  involve  model  tensors. 
Sometimes,  these  tensors  are  isotropic  and  governed  by  a  single  parameter.  The 
viscous  terms  account  for  low  Reynolds  number  effects  and  strong  inhomogeneity. 
They  do  not  involve  any  additional  parameters  and  were  derived  via  analysis  and  the 
condition  that  the  model  be  exact  as  it  approaches  a  wall  (in  the  laminar  sub  layer). 

The  zeroth  moment  of  the  PDF  equation  (Equation  (4))  is  the  mass  conservation 
equation.  The  first  velocity  moment  of  the  PDF  equation  gives  the  momentum 
equation, 


(5) 


This  implies  that  the  acceleration  is  given  by  an  =-pn+(/juln) The  viscous 

contribution  to  this  acceleration  is  necessary  only  if  the  viscosity  is  not  constant. 
Taking  the  moment  of  the  modeled  PDF  equation  with  respect  to  gives  the 

Reynolds  stress  transport  equation, 


(6) 


where  Tnmi  =  jv'nv'nv\  fdv  and  K  =  \Rn  is  the  turbulent  kinetic  energy.  The  tensors  Gljt 
Htj ,  and  Jtj  determine  the  model.  Complex  dissipation  and  pressure-strain  models 
can  be  implemented  via  these  tensors. 

The  equation  for  the  total  resolved  (or  mean)  kinetic  energy, 


The  resolved  kinetic  energy  correctly  loses  energy  as  a  result  of  large  scale 
dissipation,  and  via  turbulence  production.  It  is  completely  specified  and  does  not 
depend  on  the  model  coefficients.  The  details  of  these  derivations  can  be  found  in 
Chartrand  (2005). 

When  implementing  the  Fokker-Planck  collision  model  (Eqn.  4)  on  a  coarse 
mesh,  it  is  attractive  to  make  the  change  of  variables  f  =  ln(/) .  If  f  is  close  to 
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Gaussian  (which  is  expected)  then  /  will  be  close  to  parabolic.  This  parabola  can  be 
accurately  resolved  and  interpolated  by  the  three  points  available  in  our  scheme.  The 
evolution  equation  for  f  is, 


df  df  ,  v  df  _  •  .  df  d 

—  +  vi-J—Jr(ai-amA^—  =  -Gii-  G„v,^-  +  — 
dt  'etc.  V  mesh'dvi  "  ,}  1  dv.  dv. 
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n  n 


K  K  dv. 


(8) 


While  there  are  more  terms  to  compute  in  this  version,  the  equation  for  f  is  much 
more  accurate  to  solve  numerically.  In  addition,  low  order  methods  and  simple  (3 
point)  difference  stencils  suffice  because  /  is  expected  to  be  very  close  to  quadratic. 


The  models  for  the  tensors  Gv,  Htj ,  and  J0  require  a  time  scale  to  be 

dimensionally  correct.  For  this  reason  an  additional  transport  equation  for  the 
timescale  must  be  included  in  the  model.  The  un-oriented  eddy  model  therefore  used 
the  standard  epsilon  transport  equation  for  this  purpose  since  it  is  very  commonly  used 
in  RST  models  as  well.  The  oriented  model  obtains  the  timescale  from  the  orientation 
vector  whose  length  represents  the  inverse  of  the  eddy  size.  So  no  scale  equation  is 
necessary  in  the  most  recent  model  implementation. 

The  Original  Model 


The  original  collision  model  first  proposed  by  Perot  and  Chartrand  (2005)  was 


g„  =  c;2  s0  +  c”  Wij+±csp2—s, 


P'1 


R  R 


K 


CAj 


(9) 


jv = -f  kc;a 


(10) 

(11) 


/  K ^  is  the  modified  dissipation  that  goes  to  zero  in  regions 

of  strong  inhomogeneity  such  as  near  walls,  and  P  =  -Rnmunm  is  the  standard  turbulent 
production  rate.  The  frame  invariant  strain-rate  and  rotation-rate  tensors  are 


where  e-e 


l  +  10v 
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respectively  StJ  =~(w,.7  +uj.,)  and  =  ~(.utJ  ~  uj,,)  +  £tjSik  >  where  Q*is  the  rotation 
rate  of  a  non-inertial  frame  of  reference. 


For  comparison  with  classic  RST  models,  the  equivalent  Reynolds  stress 
transport  equation  would  be, 


— ~  ^ - W-/?  H - T  +  ( u  R  +  u  jR  .  1  — 

dt  dx,  dx.  V  ,J  J  ,J  j  ) 


+  C’p2jRm„-2-^-RmRs 


V  P 


+  4-KC'S 


8  8R 


mn 


3  p2  *" 

L/A  ;  L/A  ; 


2v 


dx. 


mn 


v  ^  y 


(12) 


Note  that  the  model  constant  Cd  does  not  affect  the  Reynolds  stress  transport 
equation.  However,  it  does  have  an  effect  on  the  higher  order  moments  (such  as  Timn) 

and  the  turbulent  transport  term.  This  constant  can  be  related  to  the  Kolmorgorov 
constant  (Pope  2000).  The  other  model  constants  are  actually  parameters  and  are 
given  by, 


Cs  = 

L>2 


v  +  v. 


IF 


/~1W  _ 

CP2  “ 


V  +  V, 


-AF  C*  =-0.2F2  +  .006— 

p2  £ 


(13) 


K2 

where  the  eddy  viscosity  is  given  by  v,  =  .12 F—r  and  F  =  ^-de\(RiJ  / k)  is  the  standard 

£ 

two-component  parameter  that  is  unity  in  isotropic  turbulence  and  zero  for  two- 
component  turbulence. 

The  transport  model  for  the  epsilon  equation  is  standard  and  is  given  by 


?f+u,^-  =  —  (CslP-Cc2e) 

dt  dx.  K  12 


d  ds 

ox,  ox. 


(14) 


where  CEl  =  1 .43 ,  Ce2  =11/6,  Ce3  =  0.83 ,  and  fairly  standard  values. 
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Inspiration  for  OEC 


The  analysis  above  lays  the  groundwork  for  the  oriented  eddy  collision  model.  A 
brief  review  of  PDF-based  models  is  in  order. 

Boltzmann  and  Fokker-Planck 

It  is  helpful  to  begin  with  a  simple  case,  and  not  consider  complications  such  as 
colliding  oriented  eddies.  Instead,  consider  a  collection  of  particles:  an  expression  can 
be  found  for  the  number  of  particles  that  have  some  velocity  v,  at  location  x,  and  time 

t ,  called  a  number  density  function.  The  more  familiar  probability  density  function  is 
simply  the  number  density  divided  by  the  total  number  of  particles  under  consideration. 
Let  /(v,,x i,t)  be  the  probability  density  function.  Using  this  function,  one  can  arrive  at 

several  useful  quantities:  multiplying  /  by  v,  and  integrating  over  all  of  the  possible 
velocities  (that  is,  taking  the  first  moment  of  /  and  integrating  over  velocity  space), 
one  can  arrive  at  the  mean  velocity  for  the  collection  of  particles,  Ui : 


(15) 


where  I  and  dvi  imply  a  triple  integral  over  v,,z  =  1,2,3.  If  the  mean  velocity  can  be 


found,  perhaps  another  quantity  of  interest,  u,Uj  can  be  found.  Taking  the  second 
moment  of  /  with  the  fluctuating  velocities  v,-Uit  recalling  the  fluctuating  velocities 
are  the  total  velocity  of  a  given  particle  v,  with  the  mean  velocity  of  all  particles  U, 
subtracted  off: 


(16) 


once  again  recognizing  that  a  triple  integral  exists  in  Equation  (15).  Equations  (15)  and 
(16)  represent  the  statistical  mechanics  of  the  collection  of  particles  but  say  nothing 
about  the  physics  present  in  that  /  has  yet  to  be  prescribed.  One  of  the  simplest  ways 

to  describe  the  time  evolution  of  a  PDF  is  through  the  Boltzmann  equation,  which 
essentially  models  particle  collisions  by  relaxing  their  PDFs  to  the  mean 
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(17) 


d  r  d  . 

— f(y, ,  x,. ,  t)  +  v,.  —  /(v, ,  x, . ,  1) 

Ot  OX; 

+«i  x"  7(v!  •  xi .  0 = 4-  fiy,  ,x„0  L 

ov(.  3/ 

with  a,  representing  some  body  (external)  force  that  may  be  acting  on  the  fluid  (such 

as  a  Coriolis  term)  and  the  right  hand  side  representing  the  way  in  which  the  average 
of  all  collisions  over  time  affects  the  PDF.  The  left  hand  side  of  Equation  (17)  is  exact, 
while  the  right  hand  side  is  that  which  requires  a  model,  the  so-called  "collision"  term. 
An  approximation  to  the  Boltzmann  equation,  such  as  the  one  originally  proposed  by 
Bhatnagar,  et.  al  (1954)  can  be  employed  and  slowly  brings  /  to  an  equilibrium  value, 
which  usually  means  a  Gaussian  distribution.  Once  a  form  of  /  has  been  chosen,  it 
can  be  used  in  Equation  (15)  and  the  mean  velocity  found  (the  method  in  which  this  is 
done  will  be  discussed  later).  A  linear  relaxation  model  may  also  be  employed  (see 
Perot  &  Chartrand’s  earlier  work).  Interestingly,  for  low  density  flows  (meaning  flows  in 
which  few  particle  collisions  occur),  a  simple  collision  model  returns  the  ideal  gas  law, 
the  viscous  terms  of  the  Navier  Stokes  equations,  Fourier  heat  conduction  and  many 
other  physical  processes.  Thus,  this  method  is  suited  for  Newtonian  flows,  but  might 
not  work  well  as  a  turbulence  model.  Inspecting  the  second  moment  and  plugging  the 
Boltzmann  equation  in  to  Equation  (16)  yields  an  unfortunate  result:  this  simple 

collision  model  predicts  that  the  Reynolds  stresses  are  zero,  duiuJ/dt  =  0.  Not  only  is 

this  approach  flawed,  it  is  in  fact  useless  for  capturing  the  behavior  of  a  turbulent  flow. 
This  is  due  to  the  fact  that  the  Boltzmann  equations  look  at  fluid  interactions  purely  as 
a  viscous  phenomenon  with  a  single  relevant  time  scale,  an  idea  which  sounds  familiar 
from  previous  turbulence  models  considered.  It  was  already  determined  that  this 
assumption  will  never  capture  turbulence  properly,  and  it  is  no  surprise  that  this 
method  fails. 

An  alternative  to  the  Boltzmann  equation  is  the  Fokker-Planck  (FP)  equation 
(also  referred  to  as  the  Kolmogorov  forward  equation  from  Kolmogorov  (1942)),  which 
describes  the  time  evolution  of  a  PDF  in  a  more  complicated  way: 

5  5  5 

— f(vi,xi,t)  +  vi—  f(v.,xi,t)  +  ai—f(vi,xi,t) 
ot  oxi  ovi 

d  d2  (18) 

=  -a  —  (v. Ut  )f(vi ,  v,. ,  0  +  f3  — T  /  ( v,. ,  x . ,  t) 

ov.  dv. 

with  a  and  fi  model  constants.  Equation  (18)  (adapted  from  Perot  &  Chartrand  2005) 
represents  one  of  the  simplest  Fokker-Planck  collision  models  with  the  right  hand  side 
of  Equation  (17)  replaced  by  two  terms.  The  right  hand  side  of  Equation  (18)  must  be 
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generalized  in  order  to  be  employed  for  PDF  turbulence  model  methods.  Pope  and 
others  proposed  a  generalized  form,  and  use  it  extensively  in  their  PDF  method  work 
(see  Pope  2000,  Pope  1982,  Pope  1983,  Perot  &  Chartrand  2005  and  more): 


3  3  3 

— /(v„*„0  +  v,—  f{v„x„t)  +  a—  f(v„X',t) 
dt  ox ,  <7v, 


(19) 


where  Gy  is  a  tensorial  modeling  parameter,  v  the  fluid  viscosity,  and  noting  the 

addition  of  a  second  order  spatial  derivative  (Laplacian)  of  the  PDF.  Unlike  the 
Boltzmann  equation,  the  FP  approach  captures  more  than  just  viscous  Navier-Stokes, 
and  the  second  moment  is  not  zero. 

The  two  methods  mentioned  above  relate  to  the  aforementioned  Langevin 
equation:  a  Langevin  approach  can  be  employed  to  solve  the  resulting  PDF  transport 
equations  arrived  at  from  plugging  either  the  Boltzmann  or  the  Fokker-Planck  collision 
models  in  to  Equation  (17)  and  the  resulting  expression  for  /  in  to  Equations  (15)  or 
(16).  This  is  because  the  equations  may  be  solved  using  "normal"  methods,  that  is 
using  a  finite-difference  or  finite-element  method  or  may  be  solved  by  using  a  particle 
approach,  the  details  of  which  will  be  avoided  here.  Using  a  Langevin  approach  makes 
no  changes  to  the  underlying  physics  -  it  is  simply  a  particle  method  solution.  When  a 
Langevin  method  is  used  with  a  Boltzmann  equation,  this  is  often  referred  to  as  a 
Lattice-Boltzmann  method  (as  the  Boltzmann  equation  is  solved  in  a  lattice  of  points, 
dictated  by  the  Langevin  equation).  Solving  a  PDF  collision  method  in  this  way  is  akin 
to  solving  Navier  Stokes  without  turbulent  terms.  Solving  the  Fokker-Planck  model 
using  a  Langevin  method  results  in  a  means  of  solving  the  Reynolds-averaged  Navier 
Stokes  equations  in  the  form  of  a  Reynolds  stress  transport  model.  This  method  is 
referred  to  by  Pope  as  the  Generalized  Langevin  Method  (GLM)  (Pope  2000,  Pope 
1994,  Haworth  1985).  Various  forms  of  the  Fokker-Planck  model  lead  to  various  forms 
of  Reynolds  stress  transport  models,  ranging  from  the  simpler  Launder,  Reece,  and 
Rodi  (Launder,  et  a/.  1974)  to  more  complex  forms.  The  fact  that  well-known 
turbulence  models  emerge  from  the  steps  above  could  be  considered  affirmation  that 
the  analysis  was  correct.  However,  simply  returning  to  a  statistical  mechanics-based 
version  of  a  well  known  turbulence  model  family  also  means  that  the  new  PDF  method 
inherits  many  of  the  previous  problems  associated  with  RST  models,  most  important  of 
which  is  the  inability  to  capture  linear  (rapid  distortion  theory  limit)  turbulence.  Despite 

the  fact  that  PDF  methods  require  no  model  for  the  triple  correlation  ukuluJ  (assuming 

they  remain  in  PDF  and  not  RST  form),  they  still  suffer  from  many  problems.  This  calls 
in  to  question  the  need  to  accurately  model  the  triple  correlation  and  suggests  that 
perhaps  it  is  in  the  pressure  term  that  the  missing  physics  may  be  found. 


ii 


Advanced  methods 


Taking  a  step  back  from  PDF  methods  for  a  moment,  the  recent  work  of 
Reynolds  &  Kassinos  (1996)  will  be  considered  briefly.  Reynolds  &  Kassinos  wished  to 
capture  rapidly  deforming  homogenous  turbulence  with  a  Reynolds  stress  transport 
model.  They  hypothesized  that  the  stress  tensor  was  not  enough  to  capture  rapid 
distortion  theory  limit  turbulence  as  information  about  the  turbulent  structure  was  not 
contained  within  such  a  quantity.  Among  others,  they  proposed  a  general  model  which 
transports  a  single,  rank  two  tensor,  the  “eddy  axis  tensor”  which  characterizes  the 
shape  and  orientation  of  a  turbulent  eddy.  The  model  employed  algebraic  equations  of 
state  (as  opposed  to  a  stress  tensor)  and  two  scalar  quantities  thus  containing 
information  about  the  dimensionality  and  “componentality”  of  the  turbulence  (Reynolds 
&  Kassinos  1996).  The  model  managed  to  capture  many  linear  turbulence  cases 
exactly,  which  was  the  first  ever  demonstration  of  an  RST-like  turbulence  model 
providing  accurate  solutions  in  this  limit  of  turbulence.  Unfortunately,  many  considered 
the  model  difficult  to  understand  (this  author  included)  mainly  due  to  the  fact  that  the 
entire  model  was  formed  in  wave  (Fourier)  space.  Despite  this  issue,  the  work 
produced  a  powerful  idea:  Perhaps  the  failure  of  PDF-based  turbulence  models  lie  not 
in  the  formation  of  the  PDF  collision  model  ( e.g .  Fokker-Planck)  but  instead  in  a 
previously  unimagined  missing  unknown,  perhaps  related  to  orientation  of  turbulent 
eddies  or,  in  the  case  of  PDF  methods,  fluid  particles. 


Perot  &  Chartrand  (2005)  proposed  a  very  general  Fokker-Planck  collision 
model  with  more  unknowns  than  the  general  Fokker-Planck  model  proposed  by  Pope 
(2000): 
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(20) 


with  an  additional  term  added  to  the  end  of  Equation  (20)  to  account  for  mesh  motion, 
which  was  related  to  a  numerical  method  used  by  Perot  &  Chartrand  to  solve  their 
generalized  Fokker-Planck  method  using  an  adaptive  three-point  mesh  in  velocity 
space  (Perot  &  Chartrand  2005).  In  Equation  (20),  Gtj , HtJ  and  Jtj  are  tensorial 

modeling  terms,  «,  is  the  physical-space  velocity  gradient  and  Kn  the  physical-space 

gradient  of  the  kinetic  energy.  Van  Slooten  and  Pope  (1996),  and  more  recently  Perot 
and  coworkers  have  attempted  to  overcome  the  inherent  limitations  of  Fokker-Planck 
based  PDF  turbulence  models.  Work  by  Perot  determined  that  any  extension  of  the 
Fokker-Planck  model  would  simply  result  in  a  slightly  more  complex  but  still  inherently 
limited  RST  model,  unless  orientation  was  given  to  the  fluid  particles,  that  is  a  Fokker- 
Planck  collision  model  was  formed  for  something  like  rods  or  disks  (later  called  eddies) 
rather  than  particles  which  are  spheres  and  have  no  orientation.  This  can  be  achieved 
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by  adding  an  '"extra"  unknown  to  a  Fokker-Planck  like  collision  model  yielding 
derivatives  with  respect  to  time,  space,  and  the  extra  term,  which  could  be  thought  of 
as  eddy  orientation. 

Van  Slooten  and  Pope  (1996)  furthered  the  ideas  presented  by  Reynolds  & 
Kassinos  (1995)  first  by  applying  them  to  a  PDF-based  method  solved  with  a  particle- 
based  approach  (a  Monte-Carlo  solution),  and  then  using  this  new  method  to  simulate 
inhomogeneous  linear  turbulence.  They  implemented  this  extra  information  via  a  joint 
PDF  of  velocity  and  a  "wave  vector"  which  is  related  to  the  unit  wave  vector  tied  to  a 
given  turbulent  eddy  size.  The  collection  of  these  vectors  are  referred  to  as  the 
directional  spectrum.  This  was  a  major  step  forward  in  PDF-based  turbulence 
modeling,  but  the  method  is  both  difficult  to  understand  and  expensive  to  solve, 
requiring  a  large  statistical  sample  in  order  to  return  reasonable  results  from  the 
particle-based  solution.  Furthermore,  Van  Slooten  &  Pope  point  out  the  need  for 
improved  dissipation  models. 


Figure  4:  Box  A  illustrates  a  hypothetical  region  of  turbulent  fluid  as  a  classic  particle  collision  model,  like 
Fokker-Planck  or  Boltzmann.  The  particles  are  spheres  and  cannot  have  any  orientation.  Box  B  is  a 
schematic  of  the  same  flow  but  with  an  expanded  collision  model  that  treats  particles  as  rods  rather  than 
spheres,  thus  including  orientation  information.  Finally,  box  C  illustrates  disks  (eddies),  which  appear  to 
be  the  shape  necessary  in  order  to  capture  linear  turbulence 


Perot  and  Chartrand  picked  up  where  Van  Slooten  and  Pope  left  off,  believing  that  the 
key  to  linear  turbulence  was  indeed  the  extra  "information"  contained  within  the  wave 
vectors.  They  chose  to  add  this  information  as  a  second  derivative  to  the  generalized 
Fokker-Planck  Equation  (17): 


(21) 


dvf 


dx 
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noting  the  second  derivative  in  orientation  space,  here  denoted  by  the  unknown  vector 
qi .  This  term  acts  as  a  sort  of  advection  in  orientation  space,  and  ends  up  being 

responsible  for  the  decay  terms  present  in  the  RST  equation  form  of  the  eddy 
orientation  vector  evolution  equation.  Perot  and  Chartrand  interpreted  this  extra 
information  as  eddy  orientation  (similar  to  Reynolds'  and  Kassinos'  hypothesis),  but 
chose  to  transform  the  PDF  collision  model  in  to  a  RST  equation  form.  The  resulting 
model  was  like  a  classic  RST  model  but  had  extra  information  inherent  to  it,  resulting  in 
a  model  which  could  capture  fast  pressure  strain  exactly,  yielded  excellent 
experimental  agreement  in  elliptical  flows,  and  calculated  linear  turbulence  exactly.  The 
dissipation  term  still  required  model  tuning,  but  nearly  all  RST  model  issues  had  been 
resolved.  Furthermore,  unlike  the  PDF  form,  when  cast  as  an  RST  model  the  turbulent 
transport  term  required  a  model,  but  models  for  this  term  are  abundant  and  not  difficult 
to  form.  Figure  4  is  a  schematic  illustration  of  this  concept.  A  critical  difference  exists 
between  the  model  proposed  by  Reynolds  &  Kassinos  (1996)  and  that  of  Perot:  Perot's 
real-space  eddy  orientation  model  did  not  take  the  moment  of  and  subsequently 
integrate  over  orientation  space,  whereas  Reynolds  &  Kassinos  did.  This  step  allowed 
Reynolds  &  Kassinos  to  cast  their  model  in  the  form  of  a  Reynolds  stress  transport 
model  which,  in  order  to  then  incorporate  the  extra  unknown,  they  multiplied  with  the 
unknown  vector  yielding  a  third  rank  tensor  transport  equation.  By  choosing  to  forgo 
this  moment,  Perot  kept  orientation  in  the  Reynolds  stress  equation  itself.  This  enabled 
the  evolution  of  the  orientations  to  be  prescribed  in  such  as  way  that  the  full  Reynolds 
stress  transport  equation,  complete  with  included  orientation  information,  was  exact  in 
the  limit  of  linear  turbulence. 

More  recent  work  by  Perot,  Chartrand  and  Andeme  (2008)  furthered  progress 
on  the  Oriented  Eddy  collision  model,  treating  it,  for  the  most  part,  as  a  modified  RST 
model  rather  than  a  PDF  collision  model.  As  was  previously  mentioned,  Perot  & 
Chartrand  chose  not  to  integrate  over  orientation  space  thus  one  Reynolds  stress 
equation  exists  for  each  eddy  orientation  vector,  and  the  average  (that  is,  RANS-like) 
Reynolds  stress  tensor  is  a  simple  average  of  all  of  the  individual,  per-eddy  Reynolds 
stress  tensors.  Perot  and  Chartrand  proposed  the  followings  per-eddy  Reynolds  stress 
evolution  equation  (note  that  nomenclature  has  been  updated  from  Chartrand  (2005)  to 
be  consistent  with  current  versions  of  the  model): 
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with  /?<y  the  Reynolds  stress  tensor  (one  for  each  eddy),  qi  the  eddy  orientation  vector 
(and  the  quantity  that  makes  this  model  look  unusual  compared  with  a  classic  RST 

model),  turbulent  viscosity  vT  =  yjK*/Kq2  noting  an  over  bar  indicates  the  quantity 

averaged  over  all  eddies,  time  scale  1  lrR  =  -^Rq2 ,  Kronecker  delta  8y,  average 
velocity  un  and  constants  C,,C2  and  a.  Note  that  model  also  includes  information 
about  the  rotation  vector  for  a  non-inertial  frame  Qk  in  u'k ,  namely  u'k  =ul  t+ellyQk, 
with  eily  being  the  permutation  tensor.  Equation  (22)  is  the  material  derivative  of  the 

Reynolds  stress  tensor.  Equation  (23)  handles  the  stress  tensor  production,  while 
Equation  (24)  accounts  for  viscous  dissipation.  Equation  (25)  provides  a  return-to- 
isotropy  model  for  the  Reynolds  stresses  (Perot  &  de  Bruyn  Kops  2006),  and  Equation 
(26)  ensures  that  the  Reynolds  stress  tensor  and  eddy  orientation  vector  (qt )  remain 
orthogonal,  which  is  akin  to  enforcing  incompressibility  (Chartrand  2005).  The  terms 
A,  +  B ,  represent  the  return-to-isotropy  model  for  the  eddy  vectors  and  a  system 

rotation  term,  respectively.  Equation  (27)  is  of  course  the  viscous  diffusion.  It  is 
instructive  to  examine  the  eddy  orientation  evolution  equation  before  these  other 
quantities  are  described: 


-<llcUk,i 

if  —  l\ 
--  avq2+—  q, 

3  V  tr  J 

-(4+*i) 

],* 


(28) 

(29) 

(30) 

(31) 

(32) 


where  Equation  (28)  is  the  material  derivative  of  the  eddy  orientation  vector,  Equation 
(29)  handles  production  with  Equation  (30)  providing  for  dissipation  and  Equation  (31) 
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being  the  standard  diffusion  term.  Equation  (31)  has  the  eddy  vector  return-to-isotropy 
term  A,  and  rotation  term  B, . 


1 


1 +  C2  — 


[3JV„  Su]gt 


t  y 


(33) 


with  constant  C3  and  Nkj=qjqk/  q2 .  The  rotation  term  B,  ensures  that  the  model 
responds  properly  to  rotation  of  non-inertial  frame,  and  is  defined  as 


B.  =  — 


(g,/Q;)  V 

20?^  + 0.25  (q;  )! 


(34) 


with  the  vorticity  vector  defined  as  Q’k  =  eIJkuk  +  £2,  recalling  Q,  is  the  rotation  vector  for 
a  non-inertial  frame. 


The  origins  of  the  terms  above  are  described  in  detail  by  Chartrand  (2005)  They 
are  explained  briefly  below.  The  dissipation  terms  (Equations  (24)  and  (30))  come  from 
observations  of  isotropic  decay.  The  part  of  the  parenthetical  term  avq2  originates 
from  low  Reynolds  number  decay.  The  second  term  \/tr  handles  high  Reynolds 

number  decay.  These  are  arguably  the  simplest  terms  in  the  model,  and  were 
constructed  first.  Next,  the  production  terms  (Equations  (23)  &  (29))  were  constructed 
using  the  exact  linearized  Navier  Stokes  equations  for  rapid  distortion  theory  from 
Pope  (2000): 


55, 

dt 


77  ,  f  k.ki 


ui,k  + 


\ 


8;, 


2ui 


UJ .*  + 


kk, 


\ 


2uu 


R 


ki 


(35) 


dk.  ,  duk 
— -  =  -k.  — - 
dt  dxt 


(36) 


noting  that  k ,  is  Pope's  wave  vector  and  that  this  is  directly  analogous  to  the  eddy 
orientation  vector  qt .  The  similarities  between  the  exact  RDT  limit  of  the  Navier  Stokes 
equations  and  the  production  terms  for  By  and  qt  are  easy  to  see.  The  production  and 

dissipation  terms  to  achieve  viscous  RDT.  Anisotropic  decay  was  tackled  next.  It  is 
interesting  to  note  that  this  is  first  time  that  inter-eddy  interactions  must  be  accounted 
for.  Rotta's  linear  return  to  isotropy  model  (Rotta  1951)  for  the  stress  tensor  was  tested 
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first.  This  version  varied  slightly  from  Equation  (32)  in  that  rather  than  an  average 
kinetic  energy  K  being  employed,  the  local  K  =  \Rmm  was  employed.  The  so-called 

“global”  version  was  also  tested,  and  found  to  perform  better  (Perot  &  Chartrand  2005). 
Chartrand  and  Perot  also  tested  their  non-linear  “EG”  return  model  and  found  it 
deficient  in  this  application,  thus  choosing  the  "global”  Rotta-like  return  to  include  in 
OEC.  Isotropy  in  OEC  is  not  only  by  isotropic  Reynolds  stresses  (on  a  local,  per  eddy 
level)  but  also  uniform  distribution  of  unit  eddy  orientation  vectors  on  a  unit  sphere. 
Flows  tend  to  distort  their  distribution  and  a  method  is  needed  to  return  to  an  isotropic 
state.  Chartrand  (2005)  investigated  six  different  methods,  the  details  of  which  will  be 
avoided  here.  The  method  employed  in  Equation  (33)  calculates  the  normalized 
distribution  of  the  eddy  vectors  Nu  and  projects  the  eddy  vectors  according  to  the 
difference  between  the  normalizes  distribution  and  the  isotropic  normalized  distribution 
represented  by  the  Kronecker  delta  6h . 

The  Complete  Original  OEC  model 

The  original  OEC  model  evolves  two  quantities:  the  eddy  orientation  vectors,  q, 
and  the  Reynolds  stresses  Ry  with  the  kinetic  energy  k  calculated  from  the  Reynolds 

stresses  as  k  =  R„/2.  The  Reynolds  stresses  are  averaged  over  all  eddies  to  produce 
R0  which  is  then  used  in  the  momentum  equation.  The  original  OEC  model  was  posed 
as: 


f  V  •  )  =  />  -  (avq2  +  ^  -  V  +  Mtj  +  V  •  ( v  +  v, ) 

V  •  (ujq, )  =  -qkukti  -  c^(a  • vq2  +  £) &  -  ($*  +  s,)  +  V(v  +  v, )Vq, 


(37) 

(38) 


The  first  equation  above  evolves  the  eddy  orientation  vectors  q,  while  the  second 
handles  the  local  (per-eddy)  Reynolds  stresses  RtJ.  Return  to  isotropy  for  £,is  handled 
by  q,K ,  defined  as 


r  r  Up  ^ 


1  +r  Dn  v 


\yT  J 


(3^-4) 


9k 


(39) 


with  Nh .  =[(1 /q,2  and  constants  CqRUp  and  .  Note  that  CqRUp  has  been 

abandoned  in  new  versions  of  OEC,  discussed  below.  Return  to  isotropy  for  Rtj  is 
similarly  handled  by  Ru*  defined  as 
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(40) 


with  constants  CUpHR  and  C0”^.  MtJ combines  rotation,  eddy  vector  return  to  isotropy 


and  a  term  which  ensures  <7,  and  Rlt  are  orthogonal,  as  shown  in  below: 


(41) 


It  is  important  to  note  that  several  major  nomenclature  changes  have  occurred  since 
this  form  of  the  model  was  created.  Hats  -  now  indicate  models  based  on  a  variant  of 
eddy  vector  (to  be  introduced  later)  and  over  bars  now  indicate  averages  taken  over 
eddies.  The  eddy  vector  return  term  q,R  is  now  denoted  simply  ^4,  while  the  system 
rotation  term  s,  has  become  Br  To  be  consistent,  the  Reynolds  stress  return  to 
isotropy  term  RtfR  is  now  4;(to  avoid  confusing  superscripts),  the  per-eddy  local  kinetic 

energy  is  K  and  the  eddy-averaged  (often  called  global)  kinetic  energy  is  simply  K. 
The  same  scheme  is  used  for  the  Reynolds  stresses  and  eddy  vectors. 

Validation 

Previous  work  on  the  oriented  eddy  collision  model  has  resulted  in  the  model 
being  validated  across  a  wide  variety  of  cases,  both  simple  and  complex  turbulent 
flows.  Considering  OEC’s  implementation  in  OpenFOAM,  numerous  validation  cases 
were  re-run.  Furthermore,  several  major  modifications  have  been  made  to  the  basic 
model  including  the  addition  of  terms  to  handle  near-wall  asymptotic  behavior  of  the 
eddy  vector  and  Reynolds  stress  tensor,  the  inclusion  of  a  near-wall  eddy  rotation  term 
in  an  attempt  to  capture  non-local  wall  effects  (discussed  later),  the  removal  of  several 
modeling  constants,  and  the  derivation  of  three  additional  forms  of  the  OEC  model 
whose  aim  it  is  to  increase  the  numerical,  temporal,  and  near-wall  stability  of  the 
model.  Considering  these  changes,  and  the  necessity  to  recreate  many  benchmark 
cases  in  a  manner  useable  by  OpenFOAM,  many  previously  performed  tests  were 
once  again  performed,  including  one  of  the  simplest  cases,  isotropic,  homogenous 
decay: 
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Table  1:  Initial  conditions  from  de  Bruyn 
Kops’  DNS  of  isotropic  homogeneous  decay 


de  Bruyn  Kops  (DNS) 

£(m/s3) 

0.782 

K(m  V) 

0.087 

v(m2/s) 

1 ,49e-  5 

ReT 

655 

Figure  5  Isotropic,  homogeneous  decay  compared  to  DNS  data  from  de  Bruyn  Kops  ,  et  at.  (1998) 

Data  from  de  Bruyn  Kops,  et  at.  (1998)  was  again  employed  for  validation.  As  was 
expected,  agreement  between  OEC  and  the  DNS  data  was  excellent,  and  more 
complex  cases  could  be  considered. 


Perfecting  the  system  rotation  term 

The  system  rotation  term  (previously  sn  now  5,)  was  re-examined  after  long¬ 
term  stability  for  shear  cases  such  as  Matsumoto,  et  at.  (1991)  came  in  to  question. 
Several  rotation  term  options  were  considered  and  their  constants  determined.  The 
current  rotation  term  is  defined  in  as 
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(42) 


with  l/rfl  ={Kq2J2  a  decay  constant,  the  eddy  orientation  vector,  Q.  =  siJkukJ  +Q,  a 

modified  system  rotation  vector,  K  the  average  (not  per-eddy)  kinetic  energy,  N  the 
number  of  eddies  in  the  simulation  and  constants  c*and  c".  Note  the  updated 

nomenclature.  The  model  was  then  subjected  to  more  complex  flow  situations  such  as 
rotating  decay,  a  mixing  layer,  and  several  shear  and  strain  cases. 


Table  2:  Initial  conditions  for  Wigeland  and  Nagib  (1978). 


Wit 

geland  &  Nagib 

A 

B 

C 

e(iti2/s3) 

14.85 

14.67 

14.94 

2.96 

3.49 

3.36 

2.77 

3.36 

22.26  1 

K(m2/s2) 

0.098 

0.0975 

0.105 

0.045 

0.0462 

0.051 

0.029 

0.033 

0.096 

v(mz/s) 

1.8e-5 

1.8e-5 

1.8e-5 

1.8e-5 

1.8e-5 

1.8e-5 

1.8e-5 

1.8e-5 

1.8e-5 

ReT 

36 

36 

41 

38 

34 

43 

17 

18 

23 

Roj 

00 

7.52 

1.78 

00 

3.77 

0.82 

00 

5.09 

2.9 

O, 

0 

20 

80 

0 

20 

80 

0 

20 

80 

Wigeland  and  Nagib  (1978)  provide  data  which  test  both  rotating  and  non-rotating 
decay,  as  seen  in  Figure  6! 
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time  (s) 

Figure  6:  Rotating  and  non-rotating  decay  from  Wigeland  and  Nagib,  case  A. 


time  (s) 

Figure  7:  Rotating  and  non-rotating  decay  from  Wigeland  and  Nagib,  case  B. 


Figure  8  Rotating  and  non-rotating  decay  from  Wigeland  and  Nagib,  case  C 


As  Figure  6  through  Figure  8  show,  agreement  is  good  across  a  range  of  Rossby  and 
Reynolds  numbers.  Rotating  decay  was  also  tested  with  results  from  Jacquin,  et  al. 
(1990).  Note  that  only  the  highest  Reynolds  number  case  is  shown  here  as  agreement 
at  lower  Reynolds  numbers  was  excellent  and  tested  extensively  previously. 
Agreement  with  Jacquin’s  high  Reynolds  number  data  was  excellent. 

Table  3:  Initial  conditions  for  Jacquin's  rotating 
decay. 


Jacquin 

A 

B 

C 

z(m2l s3) 

11.73 

16.43 

30.93 

K(m*ls*) 

0.153 

0.288 

0.444 

v(m2/s) 

1  51e-5 

1 ,51e-5 

1.51e-5 

50 

CD 

H 

127 

281 

457 

Roj 

1.22 

0.91 

1.10 

22 


time  (s) 

Figure  9,  Rotating  and  non-rotating  decay  from  Jacquin,  ef  al.  with  RoT  =1.10  (case  C). 

Finally,  data  taken  from  Mansour,  Cambon,  and  Speziale  was  also  used  to  determine 
the  performance  of  OEC's  rotation  terms: 


0  0.25  0.5  0.75  1  1.25  1.5  1.75  2  2.25  2.5 

time  (s) 

Figure  10  OEC  compared  to  rotating  isotropic  decay  data  from  Mansour,  Cambon,  and  Speziale  (A). 
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Figure  11:  OEC  compared  to  rotating  isotropic  decay  data  from  Mansour,  Cambon,  and  Speziale  (B). 


time  (s) 

Figure  12:  OEC  compared  to  rotating  isotropic  decay  data  from  Mansour,  Cambon,  and  Speziale  (C). 
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Figure  13.  OEC  compared  to  rotating  isotropic  decay  data  from  Mansour,  Cambon,  and  Speziale  (D) 

Table  4:  Initial  conditions  for  Mansour,  etal  rotating 

cases. 


Mansour,  Cambon  &  Speziale 

A  B 

C  D 

E(mz/s3) 

0.93 

0.95 

K(mz/sz) 

0.964 

0.977 

v(m^/s) 

3.67e-2 

1 ,49e-2 

ReT 

27.2 

67.1 

Roy 

0.37  1  0.037 

0.24  0.1 

The  model  agreed  well  with  available  data  (mostly  from  direct  numerical  simulations) 
and  mixing  layer  data  from  Winckelmans,  Jeanmart,  and  Carati  (2002)  was  used  to  test 
the  model’s  ability  to  capture  the  decay  of  kinetic  energy  and  dissipation  which  differs 
spatially.  Kinetic  energy  results  are  shown  in  Figure  14  and  dissipation  results  in  Figure 
15  .  Carati’s  data  is  unique  in  a  sense  that  we  have  access  to  both  the  kinetic  energy 
and  the  dissipation  rate.  OEC’s  prediction  of  the  decay  of  both  the  average  kinetic 
energy  K  and  average  (calculated)  dissipation 


£ 


~N 


_ 3 

Y,(<l2K)va  +  Ki 


(43) 


is  quite  close  to  Carati’s  at  t  =  0.171  seconds.  However,  the  model  seems  to  slightly 
over  predict  the  kinetic  energy  at  the  later  time  and  under  predict  the  dissipation. 
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Figure  14. 
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OO  Carati  - 1  =  Os 
OO  t  =  0  071s 
OO  t  =  0  191s 
- OEC 


q qOOqo O 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5  5.5  6  6.5 

y 

Decay  of  kinetic  energy  versus  position  at  three  different  times  from  Carati,  et  al.  (2002). 


time  (s) 


Figure  15:  Decay  of  dissipation  versus  position  at  three  different  times  from  Carati,  et  al.  (2002). 


26 


r 


OEC  in  the  Rapid  Distortion  Theory  Limit 

The  addition  of  orientation  information  to  OEC  enables  the  model’s  unique 
ability  to  capture  turbulence  in  extreme  circumstances,  such  as  those  described  by 
rapid  distortion  theory  (RDT).  Pope  (2000)  covers  this  subject  in  detail.  Amongst  the 
RDT  cases  considered  and  used  for  validation  were  the  following:  Axisymmetric 
expansion,  akin  to  an  expansion  in  a  wind  tunnel  in  transverse  directions;  axisymmetric 
contraction  in  which  the  turbulent  flow  is  contracted  in  the  transverse  directions,  plane 
strain,  and  finally  shear.  The  four  cases  are  summarized  in  Table  5: 


Table  5:  RDT  cases  used  for  testing  OEC  in  FOAM. 


Axisymmetric 

contraction 

Axisymmetric 

expansion 

Plane 

Strain 

Shear 

*.« 

5 

-25 

5 

0 

^22 

-Is 

2 

5 

-5 

0 

*33 

-is 

2 

5 

0 

0 

^12 

0 

0 

0 

5 

5  =  (2  SySv)m 

n/3 5 

2n/35 

-25 

25 

0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2 

Sxt 

Figure  16:  OEC  subjected  to  plane  strain  and  compared  to  theory  (Pope  2000). 
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3 

2 


OO  Theoretical  -  R„  /  K° 
~  OO  R22  /  K° 

OO  R»  /  K° 

- 0  5*exp(X) 

—  OEC 
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Figure  17:  OEC  subjected  to  axisymmetric  expansion  and  compared  to  theory  (Pope  2000). 


Figure  18:  OEC  subjected  to  axisymmetric  expansion  and  compared  to  theory  (Pope  2000) 
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Figure  19:  A  closer  look  at  the  behavior  of  the  normalized  stress  component  R22  compared  to  theory 

(Pope  2000). 


OEC  in  Shear  Flow 

After  testing  each  term  in  the  model,  OEC  was  subjected  to  homogeneous 
turbulent  shear  flow  in  order  to  compare  with  data  available  from  Matsumoto,  Nagano, 
and  Tsuji  (1991)  as  well  as  L.  Le  Penven,  J.  N.  Gence,  and  G.  Comte-Bellot  (1985). 
This  provided  a  means  of  testing  the  model  as  a  whole  while  remaining  geometrically 
simple  and  not  requiring  wall  boundary  conditions. 

Table  6:  Summary  of  shear  flow  cases  used  to  validate  OEC 


Matsumoto 

LePenven  A 

LePenven  B 

SK/e 

30.6 

4.71 

0.43 

0.33 

R©t 

18.18 

152 

612 

846 

Strain 

Tensor 

'0  28.28  0^ 

0  0  0 

,0  0  0y 

'0  30  0^ 

0  0  0 

,0  0  ov 

^5.48  0  0  ' 

0  1.99  0 

,  0  0  -7.47  j 

"8.86  0  0  > 

0  -2.36  0 

,  0  0  6.50; 
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time  (s) 

Figure  20:  Comparison  to  data  from  L  Le  Penven,  J  N  Gence,  and  G.  Comte-Bellot,  case  A. 


Figure  21:  Comparison  to  data  from  L.  Le  Penven,  J  N.  Gence,  and  G  Comte-Bellot,  case  B. 

Results  for  Matsumoto,  et  al.  (1991)  low  Reynolds  number  flow,  Rer  =  18 ,  are  shown  in 
Figure  22: 
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Figure  22:  Shear  data  at  a  turbulent  Reynolds  number  Rer  =18  from  Matsumoto,  Nagano,  and  Tsuji 

who  provide  data  for  the  time  evolution  of  the  anisotropy  tensor.  Time  in  non-dimensionalized  by  the 
shear,  S. 

The  model  agrees  quite  well  with  the  data  provided  from  Matsumoto,  et  al.  (1991). 
Results  for  higher  Reynolds  number  ReT=  152  flow,  for  a  much  longer  span  of 
characteristic  time  scale  St,  is  shown  below  in  Figure  23: 


Figure  23:  Shear  data  at  a  turbulent  Reynolds  number  Rer  =152  from  Matsumoto,  Nagano,  and  Tsuji 
who  provide  data  for  the  time  evolution  of  the  anisotropy  tensor  4, . 
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Again,  agreement  is  excellent  with  Matsumoto’s  data.  With  the  successful 
benchmarking  of  OEC  with  Matsumoto’s  shear  case,  the  model  was  ready  for  further 
development.  Of  course,  solid  boundaries  are  always  a  concern  with  RSTM-like 
models,  and  are  currently  a  subject  of  intense  research.  Before  tackling  that  problem, 
however,  the  problem  of  temporal  stability  will  be  addresses. 

Temporal  Stability 

Increasing  the  temporal  stability  of  the  code  was  addressed  next  and  a  3rd  order 
low-storage  Runge-Kutta  time  marching  scheme  was  successfully  implemented  and 
tested  in  OpenFOAM.  For  temporal  discretization,  the  program  utilizes  a  three  step 
Runge-Kutta  time  marching  method  (RK3),  which  is  second  order  accurate.  Denoting 
intermediate  solution  steps  as  y  and  | ,  we  arrive  at  the  following  low  storage,  second 
order  accurate  form  of  the  hybrid  RK3  found  in  the  code: 


\ 


yn+\  =  yn +{&}/  y  ■ 

V  ”+2  J 


(44) 


where  y  ,  represents  the  intermediate  velocity  (flux),  pressure,  Reynolds  stress,  eddy 


n+— 


vector,  or  kinetic  energy  information.  The  first  step  of  RK3  uses  the  explicit  Euler 
method  to  arrive  at  a  solution  at  one  half  the  time  step.  The  code  then  uses  this 
midpoint  solution  to  leapfrog  to  the  end  of  the  time  interval.  Finally,  it  performs  another 
Euler  step  to  arrive  at  a  solution  at  the  next  time  step.  The  low  storage  method  trades 
off  accuracy  for  minimal  storage.  Only  two  arrays  need  be  stored  for  any  given 
calculation,  the  solution  from  the  previous  step,;y„,  and  the  result  of  the  previous 

intermediate  step  y  or  p.  Implementing  this  method  in  to  OpenFOAM  was  done  for 
OEC  specifically,  and  not  in  a  general  form.  OpenFOAM  allows  for  runtime  selection  of 
time  stepping  schemes,  but  is  currently  limited  to  one  step  schemes  such  as  explicit 
Euler  or  second  order  schemes  such  as  Crank-Nicholson.  As  such,  any  higher  order 
scheme  such  as  RK3  must  be  added  to  FOAM. 

As  was  the  case  with  boundary  conditions  (discussed  below),  there  are 
generally  two  approaches  to  implementing  a  new  feature  in  FOAM.  The  first  is  to  spend 
the  effort  of  creating  a  generic,  templated  entity  (in  this  case  a  time  derivative  scheme), 
fold  this  code  in  to  the  existing  framework,  and  then  call  the  method.  While  more 
attractive  to  the  general  user,  the  time  required  to  do  this  is  often  not  worth  the  reward. 
In  this  case,  the  RK3  scheme  was  “hard  coded”  in  to  OEC  and  employ’s  FOAM’s 
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explicit  Euler  time  derivative  scheme  for  temporal  sub-stepping.  FOAM  stores  the 
previous  values  for  a  given  entity  (such  as  an  eddy  vector)  making  implementation 
easier.  Old  values  can  be  easily  recalled,  and  in  certain  circumstances  FOAM’s  default 
behavior  can  be  overridden  using  the  u.storeOldTime()”  function.  This  was  especially 
useful  when  constructing  the  RK3  scheme  in  OpenFOAM  as  there  are  a  good  number 
of  intermediate  arrays  to  be  stored  for  the  Reynolds  stress  tensor,  eddy  vector,  kinetic 
energy,  and  velocity  at  each  cell  for  every  eddy.  Future  work  on  OEC  in  OpenFOAM 
may  include  the  development  of  a  general  Runge-Kutta  time  marching  scheme  for 
users  of  the  model  and  the  general  public. 

Solid  Boundaries 

After  temporal  stability  issues  were  overcome,  effort  shifted  to  slip  and  no-slip 
boundary  conditions.  These  are  imperative  for  wall-bounded  flows  or  flows  over  objects 
which  are  of  paramount  importance  to  engineers  and  are  the  focus  of  our  current  work. 
The  model  is  currently  being  tested  using  benchmark  cases  such  as  turbulent  Couette 
flow,  turbulent  channel  (Poiseuille)  flow  and  a  backward  facing  step  In  addition,  more 
complex  cases  such  as  flow  around  an  oblate  spheroid  have  been  considered. 
Preliminary  work  on  turbulent  flow  (/?e  =  5000)  over  an  oblate  spheroid  using  OEC 
running  in  parallel  on  four  processors  is  shown  in  Figure  24  and  Figure  25  : 
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Figure  24:  Velocity  streamlines  from  turbulent  flow  over  an  oblate  spheroid  using  the  OEC  model.  The 
spheroid's  surface  is  no-slip  while  the  domain  walls  are  zero  gradient. 
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Figure  25:  Velocity  streamlines  from  turbulent  flow  over  an  oblate  spheroid  using  the  OEC  model.  The 

wake  can  be  easily  seen  in  this  slice. 


Trials  continue  with  simulations  over  solid  boundaries  including  turbulent  Couette  flow, 
turbulent  channel  flow,  and  flows  over  an  oblate  and  prolate  spheroid. 

Modifications  for  near-wall  stability 

Unfortunately,  as  is  the  case  with  many  Reynolds  stress  transport  models 
(RSTM),  problems  arise  when  walls  are  introduced  as  both  the  Reynolds  stress  and 
kinetic  energy  become  zero  at  the  boundaries.  In  order  to  maintain  stability  at  walls,  the 
OEC  model  was  recast  to  evolve  the  Reynolds  stresses  normalized  by  the  kinetic 
energy,/?’^  =RV/K .  This  necessitated  an  evolution  equation  for  the  kinetic  energy  K 

as  well.  The  eddy  orientation  vector  equation  was  unchanged.  The  new  casting  of  the 
OEC  model  became: 


?E-+V.(uJK)  =  -KutJ-'-(aW+i;)K 

+  P*K  -  A*K  +  M*K  +  V(v  +  vt  )VK 

-{< -KA')+iMl -Km-)+v-(v+v,)vr; 


(45) 


(46) 
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noting  that  the  equation  for  q ,  is  unaltered  and  that  Cp=  2.  The  first  equation  above 
evolves  the  per-eddy  kinetic  energy  AT.  In  Equations  (45)  and  (46),  P’  is  the  same  as 
Py  except  it  involves  A*  as  opposed  to  Rtj  which  is  true  for  M*  as  well.  P*  =  0.5 P*  and 
similarly  M*=0.5M*.  The  Reynolds  stress  return  to  isotropy  term  4’  is  slightly 
different  from  the  one  found  in  the  original  OEC  model,  as  shown  in  below: 


(47) 


noting  that  the  average  kinetic  energy  K  is  normalized  by  the  local  (per-eddy)  kinetic 
energy  K  Evolving  P*  allows  the  per-eddy  Reynolds  stress  P,;to  be  calculated  via 

Ry  =  R'tjK  which  does  not  present  problems  when  K  =  0.  Once  complete,  the  new^( ,  K, 
R'y  (“qkRStar”)  casting  of  OEC  was  tested  using  the  same  cases  that  were  employed 
for  the  original  qt ,  Rtj  (“qR”)  model.  Results  from  the  two  models  matched  closely. 

Near-wall  behavior  of  turbulent  eddies 

In  addition  to  ensuring  stability  at  walls,  it  is  necessary  to  ensure  that  the  eddies 
(i.e.  the  eddy  orientation  vectors  and  subsequently  the  per-eddy  Reynolds  stresses) 
interact  with  the  region  near  a  wall  (the  large-scale  damping  effect).  In  addition,  the 
eddies  must  align  themselves  properly  to  ensure  they  are  not  embedded  within  the 
wall' 


Figure  26:  Eddies  that  intersect  solid  boundaries  must  be  rotated  out  of  the  way.  A)  This  rotation 
preserves  the  magnitude  of  the  eddy,  which  does  not  affect  the  near-wall  dissipation  B)  This  scaling 
achieves  the  same  goal,  but  affects  the  near-wall  dissipation 


As  such,  an  additional  term  is  effectively  added  to  (although  not  explicitly  stated  in) 
both  the  q,  and  RtJ  or  R*  equations.  One  method  rotates  the  eddy  vector  away  from 

the  wall  while  maintaining  its  magnitude,  shown  in  Figure  26  A.  This  method  does  not 
affect  the  near-wall  dissipation  by  maintaining  the  length  scale  (eddy  vector 
magnitude).  A  second  method,  illustrated  in  Figure  26  B,  changes  the  magnitude  of 
the  eddy  vector,  which  of  course  affects  the  eddy  vector  magnitude  (by  decreasing  it) 
and  thus  the  near  wall  dissipation.  For  the  first  case,  the  angle  between  the  old  and 
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new  eddy  vectors  is  calculated  and  Rodrigues’  rotation  formula  applied  to  align  ^  or 
/?’in  with  the  new  eddy  vector  incompressibility  (i.e.  orthogonality  between  q, and  Rv  or 
R'j).  Implementation  details  and  the  algorithm  employed  to  perform  this  rotation  are 
discussed  later. 

The  presence  of  walls  in  a  turbulent  flow  imparts  so-called  non-local  effects  on 
the  flow,  specifically  affecting  turbulent  redistribution.  Durbin  (2001)  discusses  two  of 
the  most  common  methods  that  near-wall  modeling  is  achieved,  pressure  echo  and 
elliptical  relaxation.  Both  methods  seek  to  alter  turbulence  quantities  near  to  but  not  at 
a  wall  in  order  that  the  model  return  more  realistic  results.  Solid  boundaries  tend  to 
cause  regions  of  strong  inhomogeneity,  production  and  shear.  The  region  acts  to 
suppress  wall-normal  turbulence,  which  can  have  a  drastic  effect  on  the  nature  of  the 
near-wall  Reynolds  stress  tensor.  Unfortunately,  most  RSTMs  lack  a  mechanism  to 
ensure  this  behavior,  thus  special  consideration  must  be  made.  OEC  is  no  different, 
and  the  aforementioned  "near  wall  rotation"  of  the  eddy  vectors  is  this  model’s  novel 
solution  to  the  problem.  Great  care  must  be  taken  when  attempting  to  use  RST  models 
near  solid  boundaries  where  the  velocity  and  Reynolds  stresses  tend  to  zero.  At  the 
moment,  wall  functions  and  damping  are  the  most  popular  methods  employed  to 
handle  solid  boundaries.  Not  only  is  it  imperative  that  the  value  at  the  boundary  be 
prescribed,  but  the  model  must  also  behave  properly  as  it  approaches  the  wall, 
meaning  the  model's  asymptotic  behavior  must  be  considered.  If  the  fluctuating  velocity 
is  considered  to  be  a  smooth  function  of  the  distance  from  the  solid  boundary^ ,  then  it 

can  be  expanded  as  a  Taylor  series,  viz.  u,  =  Pt+q^+^y2  with  pn  qt ,  and  rt  functions 
of  the  wall-tangent  directions,  and  truncating  higher  order  terms.  If  the  velocity  at  the 
wall  is  zero,  «,(.y  =  0)  =  0,  then  p,=0  which  implies  «,  and  w3  (in  the  tangential  x  and 
z  directions,  respectively)  approach  the  boundary  like  y .  Furthermore,  if  continuity  is 
invoked,  it  is  found  that  velocity  in  the  wall  normal  direction  u2  approaches  the  wall 
like  >>2 .  Using  this  information,  the  near  wall  asymptotic  behavior  of  the  individual 
Reynolds  stress  tensor  components  can  be  assessed:  uxu{ ,  «3w3 ,  and  u,w3  will 
approach  like  y2 .  uxu2  and  u2u3  will  go  like  y3 ,  and  u2u2  like  yA  (Durbin  2001,  Pope 

2000).  It  is  not  trivial  to  ensure  this  behavior,  and  is  an  area  of  active  research  for  the 
Oriented  Eddy  Collision  model.  The  near-wall  region  creates  other  difficulties:  the 
highest  shear  rate  is  often  located  at  a  solid  wall,  and  the  normal  velocity  being  forced 
to  zero  at  the  wall  tends  to  affect  the  flow  away  from  the  wall  via  pressure  (often  called 
"wall  blocking"). 

Another  subtlety  which  arises  when  implementing  wall  boundary  conditions  for 
the  eddy  orientation  vectors  qt .  At  a  shear  free  (slip)  wall,  the  eddies  align  themselves 

to  be  tangential  to  (in  the  plane  of)  the  wall  and  the  magnitude  of  the  eddy  should 
remain  unchanged.  At  a  no-slip  wall,  however,  the  eddy  size  should  approach  zero. 
Considering  ql  has  units  of  1  /  length,  an  eddy  of  zero  size  would  correspond  to  a  q,  of 
infinity.  In  order  to  avoid  this  problem,  the  OEC  model  has  once  again  been  re-cast  to 
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evolve  the  eddy  length  itself,  L{=q  / q2 ,  which  has  units  of  length  and  thus  can  be  set 
to  zero  at  no  slip  walls.  This  will  be  covered  in  the  next  section. 

Avoiding  Troublesome  Boundary  Conditions 


Numerous  changes  have  been  made  to  both  the  OEC  model  and  its 
nomenclature,  most  of  which  centered  on  casting  the  model  in  such  a  way  that  it  would 
be  stable  at  solid  boundaries  and  numerically  tractable.  As  mentioned  previously, 
evolving  a  quantity  like  the  eddy  vector  qi  (hats  will  be  dropped  hereafter)  is 

troublesome  as  the  quantity  goes  to  infinity  at  a  solid  boundary  if  the  eddy  length  scale 
goes  to  zero.  This  is  a  problem  separate  from  near-wall  local  effect  discussed  in  the 
previous  section.  Since  its  initial  casting  and  subsequent  modifications  with  a 
normalized  Reynolds  stress  tensor,  several  changes  have  been  made  to  the 
nomenclature  as  well  as  a  symbolic  representation  of  the  near-wall  rotation  performed 
on  the  eddies.  The  original  eddy  orientation  evolution  equation,  in  its  current  form: 


with  near-wall  rotation  term  Wt.  The  current  form  of  the  return  to  isotropy  term  for  the 
eddy  vectors  A , : 


(49) 


with  the  isotropy  tensor  Nkj  =q,qk  /q2  determining  the  departure  of  the  eddy  vectors 
from  theory  original  (spherical)  isotropic  condition.  The  current  system  rotation  term 
5,  is  defined  as: 


(50) 


recalling  the  time  scale  is  defined  as 


(51) 


and  the  turbulent  viscosity  as 


(52) 
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Note  that,  for  simplicity’s  sake,  the  hats  (which  previously  indicated  a  “per-eddy" 
quantity)  have  been  dropped.  In  addition,  the  previous  rotation  constants  have  been 
hard  coded  as  C*  =20.0  and  C”  =0.25.  The  corresponding  Reynolds  stress  tensor  is 
currently  posed  as 
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with  near  wall  rotation  term  Wy  responsible  for  rotating  the  Reynolds  stress  tensor, 
using  Rodriguez’s’  rotation  formula,  to  be  aligned  with  the  rotated  near-wall  eddy.  The 
return  to  isotropy  term  for  the  Reynolds  stress  tensor,  now  written  AiJ .  takes  its  final 
form  as 
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noting  the  change  in  nomenclature  for  the  two  return  constants  and  CA°" .  It  is 

important  to  note  that  one  return  constant  was  eliminated,  and  Cf"  is  common  to  both 

the  eddy  evolution  equation  and  the  Reynolds  stress  equation.  In  fact,  all  model 
constants  remain  consistent  through  the  various  versions  presented  here.  A  term  is 
required  to  maintain  orthogonality  between  the  eddy  vector  and  the  Reynolds  stress 
tensor  (a  condition  which  is  also  enforced  by  the  near  wall  rotation  term  if  it  acts  upon  a 
given  eddy  vector  and  stress  tensor  in  the  domain).  This  term  is  written  as 
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with  A,  and  B,ihe  eddy  vector  return  to  isotropy  and  system  rotation  terms  defined 
above.  D  and  E  are  numerical  constants,  currently  set  to  D  =  2  and  E  =  0. 

As  was  mentioned  previously,  the  requirement  of  infinite  boundary  conditions  applied 
to  the  eddy  orientation  vectors  qt  led  to  the  re-casting  of  OEC  in  terms  of 

L,=qi /^2which  has  solid  boundary  conditions  of  Lt=  0.  Note  that  a  hat  -  now 
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indicates  a  model  quantity  based  on  L,  rather  than  q ,  .The  evolution  equation  for  the 
new  eddy  vector: 


-K-2-^)(4+s„)+i[(>'+i5r)i, 


1414 


i-+4- 


(56) 


with  fT„once  again  representing  the  near  wall  rotation  term,  which  will  be  discussed 

later.  Similar  to  the  original  form  of  the  model,  the  return  model  is  cast  as  such,  noting 
that  additional  tuning  using  high  Reynolds  number  shear  flow  from  Matsumoto  was 
employed  to  remove  a  tuning  constant  from  the  numerator: 


(57) 


A 

noting  the  isotropy  tensor  is  now  defined  as 
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and  the  turbulent  viscosity  is  cast  as 
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The  time  scale  is  now  written 
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The  system  rotation  term  for  the  Lt  -based  model  becomes 
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(61) 
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L,  above  is  an  approximation.  The  exact 


derivation  (or  conversion)  from  the  qt  -based  model  to  the  L, -based  model  returns  a 
dissipation  term  similar  to  that  in  the  original  casting  of  the  model,  namely 

if  ^  ' 


"'{er'hr  with  several  additional  terms  added  theZ,,  evolution  equation: 


A,.  +  §(>'+ 


v£2/,. 


vi2/., 


(62) 


These  terms  may  be  of  some  significance  near  solid  boundaries  and  are  the  subject  of 
future  work.  With  the  above  model  for  the  eddy  length  scale  L,,  a  corresponding  model 

for  the  Reynolds  stress  tensor,  now  based  on  Lt ,  can  be  constructed: 
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Note  the  similarities  between  the  version  of  the  stress  tensor  evolution  equation  based 
on  the  original  eddy  vector  ^and  its  current  form.  Now,  rather  than  denoting  local, 
“per-eddy”  quantities  with  a  hat,  we  denote  quantities  based  on  the  new  eddy  length 
scale  Lt  with  a  hat.  The  return  to  isotropy  of  the  Reynolds  stresses  based  on  I,  is 
written  as 
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and  the  corresponding  orthogonality  term 


(65) 

As  was  discussed  in  the  previous  section,  a  version  of  the  model  was  also  formulated 
for  the  evolution  of  the  Reynolds  stress  tensor  normalized  by  the  kinetic  energy, 

Ry  =  Ay /AT  where  the  star  notation  indicates  that  the  stresses  have  been  normalized. 
The  current  version  of  Rt' based  on  the  original  eddy  vector  qt : 
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with  the  return  to  isotropy  term  defined  as 
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Note  that  model  terms  based  on  the  normalized  Reynolds  stress  tensor  A*  also  carry 
with  them  an  asterisk  *.  Employing  a  normalized  stress  tensor 
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As  a  fourth  and  final  version  of  the  OEC  model,  a  version  which  combined  the  new 
eddy  vector  Lt  and  the  normalized  Reynolds  stress  tensor  Rj  was  created  in  hopes  that 
the  two  variations  would  provide  the  most  stability  near  solid  boundaries: 
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where  once  again  D  and  is  a  numerical  constant  set  to  D  =  2,  this  zeroing  final  term  in 
the  evolution  equation  and  avoiding  potential  numerical  stability  issues.  The  return  to 
isotropy  term  for  the  Reynolds  stresses  corresponding  to  the  normalizes  stress  tensor 
model  based  on  Z,, : 


.  .  CA 

A  —  __L 

"  f 

n 


vT  +  CAnv 


K~ 


LL.^ 

V  L 


K 

K 


with  the  orthogonality  term  written  as 
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As  was  the  case  with  the  previous  normalized  stress  tensor  variant  of  OEC,  an 
evolution  equation  for  the  kinetic  energy  K  is  required.  In  this  case,  this  equation  is 
constructed  using  the  eddy  length  scale  Z,( : 


The  above  four  versions  of  the  Oriented  Eddy  Collision  model  constitute  the  most 
current  version  of  the  model  under  development  in  OpenFOAM.  Considering  the 
difficult  nature  of  capturing  such  a  flow,  the  high  Reynolds  number  Matsumoto  case 
was  used  to  compare  the  performance  of  the  each  model: 
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Figure  27:  Taking  another  look  at  the  Matsumoto  ReT  =152  shear  case,  this  time  comparing  all  four 
model  variants  presented  above.  All  four  models  show  excellent  agreement 

As  seen  above  in  Figure  27,  the  three  additional  models  match  the  original  “qR"- 
“qkR*”,  which  normalizes  the  stress  tensor  by  the  kinetic  energy  and  has  a  transport 
equation  for  the  kinetic  energy,  “LR”  which  used  the  eddy  length  vector  7,  rather  than 

the  original  qt ,  and  “lkR*”which  employs  both  the  eddy  length  vector  L,  and  the 
normalized  stress  tensor  Rlj*  =  RiJ  IK . 


The  terms  preceded  by  D  and  E  are  relatively  new  to  the  model  (as  compared  to 
previous  versions)  and  warrant  discussion.  The  first  involves  the  gradient  of  a 
Reynolds  stress  tensor.  For  models  without  normalized  stress  tensor  these  terms  are 
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in  the  Reynolds  stress  equation  (for  both  models  based  on  qt  and  7,).  For  those 
involving  the  normalized  Reynolds  stress  tensor  Rtj* ,  the  term  of  interest  in  the 
Reynolds  stress  equation  is 
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The  “extra"  term  in  the  kinetic  energy  equation  is 
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The  terms  in  the  Reynolds  stress  equation  (Equations  (74)  and  (75))  and  by  extension 
that  in  the  normalized  Reynolds  stress  equation  (Equation  (76))  come  from  expanding 
the  last  term  in  the  original  OEC  formulation  (Equation  (12)  )  and  helps  ensure  the 
near-wall  asymptotic  behavior  of  the  model.  Note  that  D  is  often  chose  to  be  2,  thus 
eliminating  the  extra  term  in  the  RtJ*,  which  is  desirable  considering  it  can  cause 

numerical  difficulty  near  walls.  E  is  chosen  to  be  zero  in  an  attempt  to  ensure  that 
#2(the  average  eddy  vector  magnitude)  approaches  a  solid  boundary  like 
(2/a)/y2 where  alpha  is  a  tunable  constant,  usually  set  to  a  =15.0.  Note  that 

OpenFOAM  currently  does  not  support  tensors  above  rank  two,  and  the 
implementation  details  of  this  term  will  be  discussed  later. 

Implementing  OEC  in  OpenFOAM 

What  is  OpenFOAM? 

The  majority  of  the  initial  effort  in  this  project  focused  on  implementing  the 
oriented  eddy  collision  model  in  an  open  source  collection  of  computational  fluid 
dynamics  libraries  written  in  C++  called  OpenFOAM.  FOAM  is  unique  in  that  much  of 
the  mathematical  and  numerical  framework  required  to  perform  advanced  CFD  is 
already  in  place,  available  for  any  user  to  copy  and  modify  for  their  own  needs.  Despite 
having  a  vast  assortment  of  CFD-related  tools,  solvers,  and  utilities,  the  latest  version 
of  OpenFOAM  (version  1.7.1  from  OpenCFD  LTD)  has  few  Reynolds  stress  transport 
model  implemented.  In  fact,  it  has  only  two:  The  Launder,  Reece,  and  Rodi  (1974) 
model  and  a  variant,  the  Launder  Gibson  RSTM.  Adding  the  OEC  model  to  FOAM  was 
not  trivial,  as  no  other  model  currently  in  FOAM  must  account  for  two  to  three  transport 
equations  for  every  eddy  at  every  cell.  This  amounts  to  an  entire  collection  of  transport 
equations  that  must  be  carefully  handled  within  FOAM,  and  is  the  first  construct  of  its 
type  to  be  implemented  in  FOAM.  In  its  current  form,  OEC  employs  anywhere  from  22 
to  over  1,200  eddies  for  simulations.  The  number  of  eddies  available  to  the  code  is 
controlled  by  how  the  eddies  may  be  arranged  uniformly  on  a  unit  sphere. 
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fvm : :ddt (qiINT) 

-  (1 . 0/3. 0) *fvm: : laplacian (dEf f () ,  qiINT) 

+  (1 . 0/3 . 0) *fvm: : SuSp ( ( (alpha*nu () *qsq  +  tauR) ) ,  qiINT) 

-  fvc: :div(phi_,  qiINT) 

-  (  qiINT  6  fvc:: grad (U)  ) 

-  (  Ai  +  Bi  ) 

Figure  28:  FOAM  provides  a  vast  collection  of  operators  to  streamline  the  numerical  side  of 
implementing  a  turbulence  model  like  OEC. 

Figure  28  illustrates  the  power  of  OpenFOAM  in  that  the  software  provides  a  wide 
variety  of  useful  operators  which  eases  the  task  of  implementing  a  complex  model 
such  as  OEC.  The  entry  above  constructs  the  evolution  equation  for  qt,  and  is 

contained  within  FOAM’s  “fvVectorMatrix”  entity,  the  “fv"  indicating  “finite  volume”. 
Similar  entities  for  tensors,  “fvTensorMatrix”  and  scalars,  “fvScalarMatrix”  exist.  All 
terms  on  the  left  hand  side  of  the  equation  are  cast  implicitly,  and  are  part  of  the  matrix 
on  the  left  hand  side  of  the  system  to  be  solved  which  can  be  thought  of  as  Ax  =  6  with 

A  a  rank  two  tensor  (matrix)  which  must  be  inverted,  xthe  vector  of  unknowns,  and 

b  the  vector  of  known  on  the  right  hand  side.  Operators  such  as  “fvm::ddt"  easy  to 
identify:  “ddt”  takes  the  time  derivative  of  its  argument,  in  this  “qiINT”  which  is  the  eddy 
vector  for  the  current  eddy.  Note  that  transport  equations  such  as  this  are  constructed 
for  eddy  vectors,  Reynolds  stress  tensors,  and  in  some  cases  the  scalar  kinetic  energy 
for  every  eddy  at  every  cell  location  in  the  computation  mesh.  In  FOAM,  “fvm::”  casts 
the  operator  in  the  “finite  volume  method”,  which  essentially  places  the  operator  (and 
resulting  term)  on  the  left  hand  (implicit)  side  of  the  equation,  in  A.  For  example,  the 

Laplacian  operator  (used  for  the  viscous  diffusion  term)  is  cast  implicitly  for  stability 
purposes.  The  “SuSp”  operator  makes  a  decision  about  the  location  of  the  source  term 
(and  thus  whether  it  is  cast  explicitly  or  implicitly,  placed  in  ^  or  ;£4)  based  on  its  sign. 

Alternatively,  operators  may  be  cast  using  “fvc::”,  standing  for  “finite  volume  calculus”, 
which  is  an  explicit  casting.  This  can  be  thought  of  as  placing  the  resulting  term  in  b. 
For  example,  the  convection  term  is  handled  with  a  call  to  “fvc::div”,  which  performs  an 
explicit  divergence  operation  on  the  flux  “phi_”  and  the  eddy  vector.  The  eddy  vector 
production  term  -qkukl  employs  an  explicit  gradient  operator  (there  is  no  such  thing  as 

an  implicit  gradient  operator)  along  with  FOAM’s  inner  product,  “&”.  Finally,  explicit 
source  terms  such  as  the  return-to-isotropy  A1  and  rotation  term  Bt,  which  are 
constructed  beforehand,  can  simply  be  added  directly  to  the  equation. 

Initial  conditions  for  eddy  vectors 

A  variety  of  initial  conditions  for  the  eddy  vectors  qt  and  Z,,  are  available  for  use 

with  OEC.  These  initial  conditions  are  in  the  form  of  a  collection  of  vectors  which  are 
uniformly  distributed  on  the  unit  sphere.  These  vector  lists  were  originally  created  by 
Chartrand  (2005)  and  have  been  adapted  for  use  in  OpenFOAM.  The  number  of 
eddies  employed  in  a  given  simulation  is  akin  to  the  size  of  the  statistical  sampling 
space  given  to  the  underlying  probability  density  function  evolution  equation  underlying 
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the  model.  In  theory,  the  more  statistical  sample  space  (eddies)  given  to  the  model,  the 
better  representation  of  the  underlying  physics  is  returned.  This,  however,  comes  at  a 
cost,  one  which  is  brought  to  light  as  the  details  of  implementing  such  a  system  in 
OpenFOAM  are  considered.  Specifically,  a  system  is  required  by  which  an  arbitrary 
number  of  eddy  vectors  may  be  used  in  any  given  simulation.  Based  on  the  number  of 
eddies  (N),  each  cell  in  the  computational  domain  must  be  populated  with  N  Reynolds 
stress  tensors,  N  eddy  vectors,  and  N  transport  equations  for  each.  In  addition,  model 
variants  that  employ  the  normalized  stress  tensor  Rt’  =R„/ K  require  a  third  transport 

equation  for  the  scalar  kinetic  energy.  Two  to  three  transport  equations  for  each  eddy 
at  each  physical  location  in  the  computational  mesh  (/. e. ,  at  each  cell)  requires  precise 
accounting.  Pointer  lists  are  employed  for  this  purpose  in  FOAM.  For  some  number  of 
initial  eddy  vectors  N,  a  pointer  list  with  N  entries  is  constructed  for  the  eddy  vectors 
themselves,  for  the  corresponding  Reynolds  stress  tensors,  and  if  necessary  for  the 
scalar  kinetic  energy.  A  subtlety  arises  in  this  implementation,  which  will  be  discussed 
in  the  next  section. 

As  mentioned  above,  the  eddies  are  arranged  on  a  unit  sphere  and  are  thus  unit 
vectors.  Considering  the  magnitude  of  the  eddy  vectors  governs  the  dissipation,  these 
vectors  must  be  scaled  before  they  are  evolved  and  employed  in  the  evolution  of  the 
Reynolds  stress  tensor.  The  initial  eddy  vectors  are  scaled  by  the  positive  root  to  the 
following  quadratic  equation  (with  roots p)\ 


(78) 


with  A:0  and  Re/the  initial  kinetic  energy  and  turbulent  Reynolds  number,  respectively. 


Setting  the  initial  kinetic  energy  and  turbulent  Reynolds  number  in  turn  sets  the  initial 
desired  length  scale  (or  eddy  vector  magnitude),  which  can  be  thought  of  as  the  initial 
dissipation  present  in  the  flow  under  consideration.  Note  that  the  average  eddy 
magnitude  is  calculated  by 
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where  N  is  the  number  of  eddies  employed  in  a  given  simulation.  For  the  model 
variants  which  employ  the  alternate  length  scale  Lt  (the  “LR”  and  “LkR*”  models), 
Equation  (78)  can  be  replaced  by 
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again  noting  that  with  |Z-,|  =  [(1/jV)^](z,2)]2  similar  to  Equation  (79)  above.  The  per- 
eddy  Reynolds  stresses  are  scaled  by  the  user-set  initial  average  Reynolds  stress 
tensor  .K,”  once  the  eddy  vectors  have  been  properly  scaled: 


Again  for  the  for  the  “LR”  and  "LkR*”  variants,  the  initial  Reynolds  stresses  are  : 
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The  turbulent  Reynolds  number  Re7*  is  recalculated  once  it  is  employed  for  the  initial 
eddy  vector  and  stress  tensor  scaling: 
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And  for  the  “LR”  and  “LkR*”  model  variants: 
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The  somewhat  unusual  form  of  the  turbulent  Reynolds  number  formulations  in 
Equations  (83)  and  (84)  comes  from  the  fact  that  the  OEC  model  has  no  specific 
prescriptions  for  the  dissipations ,  thus  requiring  the  complex  denominator. 

Pointer  lists 

For  every  cell  in  a  computational  domain,  there  exists  a  collection  of  eddies  in 
that  cell.  For  every  eddy,  there  is  an  associated  eddy  vector  which  has  an  evolution 
equation,  an  associated  Reynolds  stress  tensor  with  an  evolution  equation,  and  a 
scalar  kinetic  energy  which  has  an  evolution  equation  if  the  “qkR*”  or  “IkR*”  model 
variants  are  employed.  This  concept  is  illustrated  in  Figure  29: 
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Figure  29:  Schematic  diagram  of  a  collection  of  eddies  that  may  exist  in  some  turbulent  flow.  Note  that 
each  set  of  eddies  exists  at  every  cell  in  the  computational  mesh. 

Three  pointer  lists,  of  length  N  (where  N  is  the  number  of  eddies  originally  seeded  in 
the  flow)  are  constructed.  One  is  populated  with  FOAM’S  “volVectorField"  entity,  which 
stores  a  single  vector  at  every  cell  location,  responsible  for  handling  the  eddy 
orientation  vectors.  A  second  contains  a  “volSymmTensorField’’  array,  which  stores  a 
six  component  symmetric  tensor  at  every  cell,  handling  the  Reynolds  stress  tensor. 
The  third  (when  needed)  is  a  FOAM  “volScalarField”  which,  not  surprisingly,  stores  a 
scalar  at  every  cell  location,  in  this  case  containing  the  kinetic  energy.  The 
aforementioned  subtlety  comes  when  considering  the  way  in  which  this  information  is 
accessed.  Considering  each  pointer  list  entry  is  assigned  to  a  specific  eddy,  operations 
performed  across  the  entire  computational  domain  are  performed  one  eddy  at  a  time 
because  the  pointer  lists  are  iterated  through  on  a  per-eddy  basis.  To  understand  this, 
imagine  selecting  only  the  large,  downward-pointing  eddy  in  Figure  29  at  every  cell 
location  and  manipulating  one  of  its  associated  quantities.  The  alternative  of  course  it 
to  pick  one  cell  (perhaps  the  center  cell  in  Figure  29)  and  select  every  eddy  at  that  cell, 
manipulating  some  eddy’s  associated  value  at  that  cell  alone.  This  has  advantages  and 
disadvantages.  Accounting  for  the  many,  many  tensors,  vectors,  and  scalars  in  any 
given  flow  is  trivial,  as  each  pointer  list  is  of  size  N,  each  entry  corresponding  to  the 
kinetic  energy  for  one  eddy  at  each  cell,  one  eddy  orientation  vector  at  each  cell,  or 
one  eddy’s  Reynolds  stress  tensor  at  each  cell.  This  makes  performing  averages  over 
all  eddies  as  simple  as  a  summation  over  all  pointer  list  entries  and  a  division  by  N. 
This  choice  makes  operations  that  must  be  performed  on  every  eddy  at  a  given  cell 
much  more  difficult,  however.  Such  operations  are  rare,  but  require  extensive  looping 
over  each  pointer  list  at  each  cell  location,  an  expensive  operation. 


Rij  [eddy] [cell]. xx()  « — 

-Component  (11  In 

V - y - ✓  4 

this  case) 

Pointer  list  with  an  1  Location  In  mesh 

entry  for  every  eddy 

Figure  30:  Using  variable-sized  pointer  lists  for  per-eddy  quantities  in  FOAM. 
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One  of  the  most  powerful  and  useful  features  of  OpenFOAM  is  the  ability  to 
access  and  manipulate  the  components  of  a  vector  or  tensor  field  across  an  entire 
mesh  (i.e.  across  all  cells  and  boundary  patches)  without  the  need  to  explicitly  access 
each  cell  location.  In  fact,  FOAM’s  namesake,  “ field  operation  and  manipulation”, 
betrays  the  power  of  this  ability  and  makes  the  implementation  of  such  a  complex 
model  much  simpler  in  C++.  Unfortunately,  this  feature  may  only  be  used  if  access  to 
one  eddy’s  components  across  the  entire  computational  domain  is  required,  and  not 
the  opposite,  where  the  component  of  all  eddies  at  a  single  cell  is  required.  In  Figure 
30,  the  pointer  list  addressing  is  illustrated.  If  it  is  sufficient  to  access  a  given  eddy’s 
components  (or  other  associated  entities,  such  as  “correctBoundaryConditions",  a 
function  that  updates  or  recalculates  a  field’s  boundary  values)  the  cell  addressing  may 
be  omitted  altogether,  greatly  increasing  the  efficiency  of  all  such  operations. 

Gradient  of  a  rank  two  tensor 

The  evolution  equations  for  both  the  “standard"  and  normalized  Reynolds  stress 
tensor,  for  both  includes  a  term  that  involves  the  gradient  of  the  Reynolds  stress 
tensor,  as  shown  in  Equations  (85)  and  (86).  This  is  a  rank  two  tensor,  and  its  gradient 
produces  a  rank  three  tensor.  Unfortunately,  as  of  OpenFOAM  1.6,  rank  three  tensors 
were  not  accommodated  for.  The  templating  was  in  place,  but  no  operators  could 
handle  such  an  entity,  including  the  gradient  operator.  As  such,  either  an  operator  must 
be  created  that  could  return  a  rank  three  tensor,  or  a  custom  function  written  that  could 
perform  the  calculation  required  in  the  model. 


(85) 


(86) 


The  first  choice,  extending  the  existing  gradient  operator  to  handle  any  rank  two  tensor 
would  require  immense  effort  (to  make  this  operator  sufficiently  general  and  interface 
with  the  existing  operator  templates  in  OpenFOAM)  and  thus  was  deemed  more  effort 
than  it  was  worth.  The  second  option,  writing  a  custom  function  to  perform  the  desired 
gradient  in  this  model  was  instead  completed.  Specifically,  the  function  was  created  to 
calculate  the  inner  product  of  the  stress  tensor  gradient  (a  rank  three  tensor)  and  the 
gradient  of  the  kinetic  energy  (a  rank  one  tensor)  which  results  in  a  rank  two  tensor.  A 
code  snippet  from  the  function  is  provided  in  Figure  31: 
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//  internal  cells 


forAll  (meah_.C(),  cell) 

( 

Rx [cell] . x ( )  -  Rtnp [ cel 1 ] . xx (); 
Rx[cell]  .y()  -  Rtmp[cell]  .xy ()  ; 
Rx[call].z()  -  Rtmp [call ] . xz ( ) ; 


gradKgradR[cell] . xx() 
gradKgradR[call] .yy() 
gradKgradR[ cell ] . zz ( ) 


-  (  gradK[cell] . x ( 
+  gradK [ cell ] . y ( 
+  gradK [cell] . z ( 

-  (  gradK[cell] . x ( 
+  gradK [cell ] .y ( 
+  gradK [cell] . z ( 

■  (  gradK [cell] .x ( 
+  gradK [ cel 1 ] . y ( 
+  gradK [cell] . z ( 


gradRx [call] 

.XX 

0 

gradRx [call] 

.xy 

0 

gradRx [call] 

.  xz 

0 

) ; 

gradRy [call] 

.yx 

0 

gradRy [call] 

•yy 

0 

gradRy [call ] 

•y* 

0 

) ; 

gradRz [call] 

.  zx 

0 

gradRz [call] 

•  zy 

0 

gradRz [call] 

.  zz 

0 

) ; 

Figure  31.  An  example  of  the  custom  function  written  for  calculating  the  gradient  term  from  Equations 
(85)  and  (86).  Note  that  looping  over  all  cell  locations  may  be  avoided  in  circumstances  when 

Future  work  on  OEC  in  OpenFOAM  may  include  the  creation  of  a  template,  generic 
gradient  operator  that  can  take  a  rank  two  tensor  as  an  input  and  return  a  rank  three 
tensor. 

The  PISO  loop 

Most  non-steady  state  solvers  in  OpenFOAM  use  the  so-called  pressure  implicit 
with  splitting  of  operators  (“PISO”)  loop  to  correct  the  pressure  term,  details  of  which 
can  be  found  in  Ferziger  and  Perid  (2002),  Anderson  (1995),  Rusche  (2002)  and 
Jasak  (1996).  The  basic  algorithm  is  contained  within  many  OpenFOAM  solvers  and  is 
augmented  with  equations  from  Ferziger  and  Perib  (2002),  Rusche  (2002)  and  Jasak 
(1996).  To  begin  with,  the  momentum  equation  is  constructed  keeping  in  mind  that  the 
flux  (of  the  velocity)  is  treated  explicitly  using  the  last  known  value  of  the  velocity.  The 
momentum  equation  is  then  solved  using  the  last  known  value  of  pressure  on  the  right 
hand  side.  This  results  in  a  velocity  field  that  is  not  divergence  free,  but  approximately 
satisfies  momentum.  The  velocity  at  some  node  P  is  obtained  by  solving  the  linearized 
momentum  equation  (Ferziger  and  Perib  2002): 


f  op”-'  ' 


dx 


t  J 


(87) 


and  is  expressed  as  (Ferziger  and  Peric  2002): 
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with  Q"~l  the  collection  of  source  terms  that  can  be  computed  using  velocity  at  time 
m- 1  in  terms  of  the  velocity  «,m~I  (where  the  superscript  m  in  u "  is  the  current  estimate 
to  the  solution  to  the  velocity  «”+1  at  time  h  +  1)  and  /  denoting  the  neighboring  points 
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from  the  momentum  equation  discretization.  Af‘  contains  the  coefficients  based  on  the 

old  at  neighboring  node(s)  /,  while  A) is  the  same  for  the  current  node  in  question. 

and  represent  the  current  solution  to  the  velocity  at  the  node  in  question  and 

surrounding  nodes,  respectively,  and  the  *  is  meant  to  indicate  that  current  solution  is 
not  the  final  solution  and  therefore  does  not  satisfy  the  continuity  equation. 

After  the  momentum  equation  is  set  up  and  solved,  the  PISO  loop  begins  and 
performs  a  given  number  of  corrector  loops.  First,  from  the  last  solution  of  velocity,  the 
diagonal  terms  are  extracted  from  the  matrix  and  the  reciprocal  is  calculated  and 
saved.  Then,  a  Jacobi  pass  is  taken  and  the  velocity  updated  (see  Jasak  1996  and 
Rusche  2002).  The  fluxes  are  then  calculated  accounting  for  the  divergence  of  the 
velocity  field  by  removing  the  difference  between  the  interpolated  velocity  and  the  flux. 
The  inlet  and  outlet  fluxes  are  then  adjusted  to  obey  continuity  Finally,  if  requested  by 
the  FOAM  user,  another  loop  begins  (inside  of  the  main  PISO  loop)  which  iteratively 
corrects  for  non-orthogonality.  This  step  ends  the  main  PISO  loop.  From  there,  the 
pressure  gradient  is  added  to  the  velocity  field,  noting  that  the  pressure  is  the  entire 
pressure  field,  not  a  correction.  To  summarize  the  algorithm  (adapted  from  Ferziger 
and  Perid  2002): 

•  To  begin  with,  calculate  fields  at  time  rn+I  at  some  node  i  using  values  u ”  and 
pn  as  the  initial  guess  for  uf  and  pn* 1 . 

•  Build  and  solve  the  linearized  momentum  equation  for  uf . 

•  Solve  the  pressure  correction  equation  to  obtain  the  intermediate  pressure. 

•  Correct  the  velocity  w™  and  pressure  /?msuch  that  satisfy  the  continuity 
equation. 

•  Repeat  this  procedure  updating  u ,m  and  pm  to  iteratively  improve  estimates  for 
uf  and  pn+i  until  the  corrections  become  small. 

•  Finally,  advice  in  time. 

Note  that  this  method  is  an  alternative  to  very  popular  fractional  step  methods,  which 
have  been  employed  in  previous  implementations  of  OEC. 

Implementation  details  concerning  solid  boundaries 

Eddies  that  lie  close  to  or  on  a  wall  are  rotated  such  that  they  are  not  embedded 
in  the  wall.  For  a  given  eddy  vector  qt ,  the  following  transformation  is  applied  until  the 
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eddy  has  been  rotated  far  enough  from  the  wall  making  sure  that  the  magnitude  of  q, 
remains  unchanged: 
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where  ycetl  is  the  distance  of  the  eddy  in  question  from  the  nearest  wall.  Note  once 
again  that  the  hats  have  been  dropped  for  convenience.  The  loop  terminates 
when|a,|2  <ycJ  where  a(=9,  xnyand  n,  is  the  unit  normal  vector  of  the  nearest  wall. 
The  same  procedure  works  if  I,  is  substituted  for  q, .  Once  an  eddy  is  rotated,  the 
corresponding  Reynolds  stress  tensor  must  also  be  rotated  in  order  to  ensure  ^and 
/^remain  orthogonal.  Rodrigues’  rotation  formula  is  employed: 


T0=Pu+(S'j-pv)cos<f>  +  Lu  sin  t  (92) 

with  <f>  the  angle  between  the  original  and  rotated  eddy  vector,  the  cosine  between  the 
old  and  new  vectors  is  defined  as  cos<f>  =  (qlu,d -q^/^q^Wq]),  and  the  sine  subsequently 

calculated  via  sin^  =  >/l-cos2^ ,  and  Py=a,aj.  The  skew-symmetric  tensor  LtJ  is 
defined  for  each  as 
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Finally,  employing  Equation  (92)  above,  Rtj  is  rotated  via  RtJ  =  7V/?y7^ .  The  same 

procedure  works  if  Lt  is  substituted  for  q , .  Several  ideas  are  implicit  to  this  selection  of 

a  near  wall  rotation  algorithm.  First,  and  most  importantly,  the  procedure  outlined 
above  does  not  alter  the  magnitude  of  the  eddy  vector  it  is  operating  on.  This 
essentially  decouples  the  near  wall  dissipation  from  the  rotation  operation.  Whether 
this  is  physical  or  not  is  still  the  subject  of  current  research.  It  may  be  the  case  that 
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such  near  wall  realignment  does  in  fact  affect  the  dissipation  and  thus  should  affect  the 
magnitude  of  the  eddy  vectors.  Independent  of  the  magnitude  of  the  eddy  vectors,  the 
Reynolds  stress  tensor  rotation  must  conserve  kinetic  energy  both  locally  and  globally, 
as  it  is  surmised  that  eddy  realignment  (at  least  realignment  that  does  not  affect  the 
eddy  vector  magnitude)  should  not  affect  the  local  or  global  kinetic  energy.  The  method 
used  to  rotate  the  per-eddy  stress  tensors  conserves  local  kinetic  energy  ( i.e .  the  trace 
of  the  stress  tensor  remains  constant)  and  validation  cases  have  shown  that  global 
kinetic  energy  (the  trace  of  the  average  stress  tensor)  is  also  preserved. 


For  the  most  part,  eddies  have  at  least  two  non-zero  components  and  theory 
associated  orientation  vectors  are  neither  perfectly  normal  to  nor  tangent  to  a  solid 
boundary,  as  illustrated  in  Figure  32  A.  Figure  32  B  illustrates  two  troublesome 
situations.  On  the  left,  the  eddy  is  perfectly  tangential  to  the  wall  and  on  the  right  the 
eddy  is  perfectly  normal  to  the  wall.  The  above  rotation  algorithm  fails  in  these 
situations.  If  the  eddy  is  normal,  the  eddy  should  never  be  rotated  as  it  cannot  be 


A) 

Y 


B) 

Y 


wall  wall 

Figure  32:  A)  Most  eddies  are  not  aligned  tangential  or  normal  to  a  solid  boundary.  B)  Some  eddies  lay 
be  normal  to  (left)  or  tangential  to  (right)  a  wall. 


wall 


intersecting  the  wall.  This  is  case  rarely,  but  must  be  accommodated.  In  the  event  of  a 
perfectly  parallel  eddy,  the  problem  becomes  more  serious.  In  this  case,  the  normal 
component  of  the  eddy  vector  (in  this  case,  the  Y  components)  is  identically  (or  close 
to)  zero.  Even  if  the  eddy  is  embedded  in  a  wall  (see  Figure  26),  the  algorithm  above 
will  fail.  Either  it  will  fail  to  rotate  the  eddy  (as  no  changed  to  the  tangential  X  or  Z 
components  can  possibly  rotate  the  eddy  away  from  the  wall)  or  it  will  spin  the  eddy 
about  its  wall-normal  (Y)  axis  forever.  Neither  case  is  desirable,  as  the  eddy  must  be 
rotated  out  of  the  wall  but  the  standard  approach  will  not  work.  Several  possible 
solutions  exist,  including  “nudging”  the  eddy  away  from  the  wall  by  forcing  the  wall- 
normal  (Y)  component  to  be  non-zero.  Of  course,  the  sign  of  the  arbitrary  non-zero  Y 
component  will  dictate  whether  the  resulting  eddy  vector  points  toward  or  away  from 
the  wall.  There  is  no  correct  answer  to  this  question,  and  as  of  now  the  OEC  near-wall 
algorithm  chooses  to  point  all  perfectly-tangential  eddies  slightly  away  from  the  solid 
boundary  they  are  embedded  in  to. 


Moving  from  the  region  near  a  solid  boundary  to  the  wall  itself,  the  boundary 
conditions  imposed  on  the  eddy  vectors  and  Reynolds  stress  tensor  will  be  discussed. 
Two  types  of  solid  boundaries  will  be  considered:  first,  a  classic  “slip”  wall  where 
surface-normal  velocities  are  forced  to  zero  (a  no-penetration  condition)  but  tangential 
components  are  undamped,  that  is  a  zero  gradient  condition  is  set  as  the  boundary 
condition.  The  appropriate  boundary  conditions  for  velocity  are  obvious,  as  are  those 
for  other  quantities  such  as  kinetic  energy,  pressure,  and  so  on.  Appropriate  boundary 
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conditions  for  the  eddy  vectors  and  Reynolds  stress  tensors  are  less  obvious,  however. 
Equation  (94)  proposes  slip  wall  boundary  conditions  for  the  Reynolds  stress  tensor  Rtj 


R„ 


I  slip -wall 


8R 


dy 


U-  =  0  Rn=  0 


dR 


21  —  i 


dy 


R22  =  0  =  0 


23 


dR 


33 


dy 


=  0 


(94) 


And  for  the  original  eddy  vectors  qt : 


I  wall 


<1\  =0 
dy 

.  <h  =  0  . 


(95) 


In  Equation  (94),  all  the  components  of  the  stress  tensor  which  involve  a  vertical  (2) 
component  are  set  to  zero,  while  those  independent  of  the  vertical  component  are  set 
to  zero  gradient  in  the  vertical  (in  this  case,  Y)  direction.  A  different  idea  is  applied  to 
the  eddy  vectors  in  Equation  (95).  At  a  slip  wall  (indeed  at  any  solid  boundary),  only  the 
vertical  component  of  the  eddy  vector  is  allowed  grow  or  shrink  (once  again  assuming 
Y  is  the  wall-normal  direction),  and  the  two  tangential  components  of  the  eddy  vector 
are  forced  to  zero.  Selecting  tensor  components  to  be  no-slip  or  zero  gradient  in  a 
certain  direction  is  somewhat  nebulous  in  OpenFOAM.  A  boundary  condition  does 
exist  which  allows  certain  components  of  a  vector  or  tensor  to  have  zero  gradient 
boundary  conditions  applied  while  others  can  have  a  fixed  value  (i.e.  zero)  condition 
applied.  The  current  implementation  of  this  boundary  condition  does  not,  however, 
allow  for  a  zero  gradient  boundary  condition  in  a  certain  direction  to  be  applied  -  the 
zero  gradient  is  applied  to  all  directions  of  a  given  component. 

The  second  case  considered  is  somewhat  simpler,  that  is  appropriate  boundary 
conditions  at  a  no-slip  wall  (which  is  often  of  most  interest  to  engineers)  where  most 
quantities  are  forced  to  zero.  Again,  many  boundary  conditions  are  obvious:  all  three 
velocity  components  are  forced  to  zero,  the  kinetic  energy  is  forced  to  zero  as  are  all 
six  components  of  the  Reynolds  stress  tensor.  The  eddy  vector  qt  is  left  in  the  form  of 

Equation  (95)  above,  whereas  all  components  of  the  eddy  vector  L,  are  set  to  zero  at 

the  solid  boundary,  essentially  forcing  all  eddies  at  the  wall  to  be  of  zero  size,  which  is 
intuitive  physically.  The  prescription  of  such  boundary  conditions,  while  convenient  in 
the  case  of  a  no-slip  wall,  causes  several  numerical  issues  that  must  be  taken  in  to 
consideration.  Specifically,  forcing  the  eddy  vectors  L,  to  be  zero  at  solid  boundaries 
can  lead  to  unexpected  behavior  (in  the  form  of  divide-by-zeros)  unless  care  is  taken 
when  implementing  terms  that  include  I, or  1}  in  the  denominator. 
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Accomplishments 


The  Oriented  Eddy  Collision  model  has  advanced  significantly  from  its  original  form. 
The  model  has  been  tuned  and  tested  with  a  variety  of  fundamental  turbulent  flows. 
The  model  has  been  parallelized,  and  implemented  in  OpenFOAM.  Furthermore, 
extensive  research  in  to  the  model’s  behavior  near  solid  boundaries  has  been 
performed.  To  summarize: 

•  OEC  has  been  implemented  in  an  unstructured,  fully  three-dimensional,  parallel 
Navier-Stokes  code  (OpenFOAM). 

•  OEC  has  been  validated  and  tested  in  FOAM  using  previously  employed 
benchmarks  as  well  as  several  new  test  cases. 

•  Feasible  boundary  conditions  have  been  developed  for  use  on  slip  and  no-slip 
boundaries. 

•  Insight  has  been  gained  in  to  the  physical  behavior  of  turbulent  eddies  near  solid 
boundaries,  and  ways  in  which  other  Reynolds-stress  transport  models  may  be 
handled  in  such  physical  situations. 

•  The  number  of  model  constants  present  in  the  code  has  been  reduced. 

•  OEC  is  in  a  position  to  be  used  by  a  wide  audience  for  a  variety  of  both 
fundamental  and  real-world  turbulent  flows. 


Future  Work 

The  Oriented  Eddy  Collision  model  is  still  under  development.  Specifically,  the 
behavior  of  the  model  near  solid  boundaries  is  still  being  investigated  and  perfected. 
Future  work  may  include: 

•  Using  OEC  in  OpenFOAM  to  determine  the  model’s  ability  to  predict  turbulent 
transition. 

•  Continuing  to  test  and  perfect  the  model’s  behavior  near  solid  boundaries. 

•  Extending  the  FOAM  cases  to  flows  of  interest  to  engineers. 

•  Optimizing  OEC’s  performance  in  OpenFOAM. 
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